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We consider the decoherence of a single localized electron spin due to its coupling to the lattice 
nuclear spin bath in a semiconductor quantum computer architecture. In the presence of an external 
magnetic field and at low temperatures, the dominant decoherence mechanism is the spectral dif- 
fusion of the electron spin resonance frequency due to the temporally fluctuating random magnetic 
field associated with the dipolar interaction induced flip-flops of nuclear spin pairs. The electron spin 
dephasing due to this random magnetic field depends intricately on the quantum dynamics of the 
nuclear spin bath, making the coupled decoherence problem difficult to solve. We provide a formally 
exact solution of this non-Markovian quantum decoherence problem which numerically calculates 
accurate spin decoherence at short times, which is of particular relevance in solid-state spin quantum 
computer architectures. A quantum cluster expansion method is developed, motivated, and tested 
for the problem of localized electron spin decoherence due to dipolar fluctuations of lattice nuclear 
spins. The method is presented with enough generality for possible application to other types of 
spin decoherence problems. We present numerical results which are in quantitative agreement with 
electron spin echo measurements in phosphorus doped silicon. We also present spin echo decay 
results for quantum dots in GaAs which differ qualitatively from that of the phosphorus doped 
silicon system. Our theoretical results provide the ultimate limit on the spin coherence (at least, as 
characterized by Hahn spin echo measurements) of localized electrons in semiconductors in the low 
temperature and the moderate to high magnetic field regime of interest in scalable semiconductor 
quantum computer architectures. 
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I. INTRODUCTION 

Quantum computation considerations have generated 
a great deal of recent interest in the \ ( j± i 2 i 2 i ££& problem 
of electron spin coherence in semiconductors. In partic- 
ular, a localized electron spin can act as a qubit (i.e., a 
quantum dynamical two-level system) for quantum infor- 
mation processing in scalable solid-state quantum com- 
puter architectures But using such spin qubits for 
quantum information processing purposes necessarily re- 
quires long spin coherence time, and therefore under- 
standing electron spin decoherence in the solid-state en- 
vironment becomes a key issue. In this paper, we de- 
velop a quantum theory for, what we consider to be, the 
most important electron spin decoherence mechanism in 
spin-qubit based semiconductor quantum computer ar- 
chitectures. The spin decoherence mechanism we con- 
sider here is the so-called spectral diffusion mechanism, 
which has a long historyjLSi&i^iii and has been much- 
studied recently2*i2*ii*i2ii^ in the context of spin qubit 
decoherence. 

To provide a physical background for the theory to 
be presented in this paper we start by considering a lo- 
calized electron in a solid, for example, a donor-bound 



electron in a semiconductor as in the doped Si:P system. 
Such a Si:P system is the basis of the Kane quantum 
computer architecture i The electron spin could deco- 
here through a number of mechanisms. In particular, 
spin relaxation would occur via phonon or impurity scat- 
tering in the presence of spin-orbit coupling, but these 
relaxation processes are strongly suppressed in localized 
systems and can be arbitrarily reduced by lowering the 
temperature. In the dilute doping regime of interest 
in quantum computation, where the localized electron 
spins are well-separated spatially, direct magnetic dipo- 
lar interaction between the electrons themselves is not an 
important dephasing mechanism Interaction between 
the electron spin and the nuclear spin bath is therefore 
the important decoherence mechanism at low tempera- 
tures and for localized electron spins. Now we restrict 
ourselves to a situation in the presence of an external 
magnetic field (which is the situation of interest to us 
in this paper) and consider the spin decoherence chan- 
nels for the localized electron spin interacting with the 
lattice nuclear spin bath. Since the gyromagnetic ratios 
(and hence the Zeeman energies) for the electron spin 
and the nuclear spins are typically a factor of 2000 dif- 
ferent (the electron Zeeman energy being larger) , hyper- 
fine coupling induced direct spin-flip transitions between 
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electron and nuclear spins would be impossible at low 
temperature since phonons would be required for energy 
conservation. 15 This leaves the indirect spectral diffusion 
mechanism-^ as the most effective electron spin decoher- 
ence mechanism at low temperatures and finite magnetic 
fields. The spectral diffusion process is associated with 
the dephasing of the electron spin resonance due to the 
temporally fluctuating nuclear magnetic field at the lo- 
calized electron site. These temporal fluctuations cause 
the electron spin resonance frequency to diffuse in the 
frequency space, hence the name spectral diffusion. The 
specific electron spin spectral diffusion process being con- 
sidered in this work is that arising from the magnetic 
dipolar interaction induced flip-flop mechanism in nu- 
clear spins. Other local flip-flop mechanisms between nu- 
clei, such as indirect exchange interactions22i2ii22 (which 
we will briefly address in this paper in the context of 
GaAs), may be easily included in our formalism. 

Spectral diffusion (SD) is, in principle, not a limit- 
ing decoherence process for silicon or germanium based 
quantum computer architectures because these can, in 
principle, be fabricated free of nuclear spins using iso- 
topic purification. Unfortunately this is not true for the 
important class of materials based on III-V compounds, 
where SD has been shown to play a major role^i^ There 
is as yet no direct (e.g., GaAs quantum dots) experimen- 
tal measurement of localized spin dephasing in III-V ma- 
terials, but such experimental results are anticipated in 
the near future. Indirect spin echo measurements based 
on singlet-triplet transitions in coupled GaAs quantum 
dot systems^ give T2 times consistent with our theoret- 
ical results. 

SD is a theoretically challenging problem because the 
temporally fluctuating random field causing the electron 
spin dephasing is a non-Markovian stochastic dynamical 
variable arising from the complex dipolar quantum dy- 
namics of a large number of interacting nuclear spins. 
Although SD is an intrinsic dephasing mechanism con- 
tributing to the so-called T2-type (i.e. "transverse") spin 
decoherence in the phenomenological Bloch equation lan- 
guage, the effect of SD on the electron spin dynamics (for 
example, on the electron spin resonance measurement) 
cannot be characterized by a simple relaxation (or de- 
phasing) time T2 associated with a "trivial" e _4 / T2 type 
of temporal loss of coherence. The non-Markovian nu- 
clear spin bath dynamics produces a rather complex elec- 
tron spin dephasing which cannot be characterized by a 
single decoherence time (except as a very crude approxi- 
mation). It is therefore more appropriate to consider SD 
in the context of a specific experimental situation, and it 
has been traditional, going back to the seminal pioneer- 
ing work of Hahn, to study SD in the context of the spin 
echo decay in pulsed spin resonance measurements iSa In 
spin echo measurements, the inhomogeneous broadening 
effects are suppressed by design, and any amplitude de- 
cay in the consecutive pulses arises entirely from the T2- 
dephasing time associated with SD (and, of course, other 
applicable T2-type transverse relaxation processes). 



It was realized a long time agoi^ that SD due to the 
dipolar fluctuations of nuclear spins often dominates the 
coherence decay in electron spin echo experiments. All 
available theories to date are based on classical stochas- 
tic modeling of the nuclear field, a Markovian theoretical 
framework which is inevitably phenomenological since it 
requires an arbitrary choice for the spectrum of nuclear 
fluctuations. Such a classical Markovian modeling is ar- 
guably incompatible with the strict requirements of spin 
coherence and control in a quantum information device. 
In addition, recent rapid experimental progress in single 
spin measurements fSii which in the near future promises 
sensitive measurements of quantum effects in spin res- 
onance, also warrants a quantum theory of SD. In the 
current work, we develop a quantum theory of spin echo 
decay in electron spin resonance experiments due to the 
SD caused by the dipolar interaction induced flip-flop nu- 
clear spin dynamics. We emphasize that, in contrast to 
all the earlier theories of SD in the literature, our theory 
is fully quantum mechanical and incorporates the dipo- 
lar nuclear spin flip-flop dynamics microscopically with- 
out making any phenomenological statistical (Markovian 
or otherwise) approximations. To the best of our knowl- 
edge, ours is the first fully quantum theory for electron 
spin SD in solid state systems. A proper description of 
coupled spin dynamics is rather difficult due to the ab- 
sence of the usual Wick's theorem for spin degrees of 
freedom. In that regard, variations of our method may 
prove rather useful, since environmental spin baths are 
ubiquitous in any device exploiting the coherent proper- 
ties of quantum spin systems. 

Our theory produces an accurate quantitative and 
qualitative prediction of the Hahn echo decay. It was 
pointed out a long time ago^&ii that the observed time 
dependence of these echoes are well fitted to the ex- 
pression exp (— i 2 ), a behavior which can be derived 
phenomenologically^ by assuming Lorentzian Brownian 
motion for the electron spin Zeeman frequency. In our 
method this approximate behavior arises naturally from 
the collective quantum evolution of the dipolar coupled 
nuclei, without any phenomenological assumption on the 
dynamics of the environment responsible for decoherence. 
Our theory reveals that the inclusion of quantum correc- 
tions to nuclear spin fluctuation increases the degree of 
decoherence at lower to intermediate time scales, as is 
best evidenced from our explanation of the existing fac- 
tor of 3 discrepancy between the Markovian stochastic 
theory^ and experimental dataS*i2*ii of spin echo decay 
in phosphorus doped silicon. 

The rest of this paper is organized as follows. In Sec.lTTl 
we provide a brief background for the spectral diffusion 
process; in Sec. IHII we formulate the theoretical problem 
of spectral diffusion in the context of Hahn spin echo 
decay; in Sec. II VI we develop and describe the quantum 
cluster expansion which is the important theoretical re- 
sult of our work; in Sec.0we describe dual perturbation 
theories in terms of a r-expansion and a dipolar pertur- 
bation expansion which together explain the convergence 



3 



Spectral diffusion of a Si:P spin 




FIG. 1: The electron of a P donor in Si experiences spectral 
diffusion due to the spin dynamics of the eveloped bath of Si 
nuclei. Of the naturally occuring isotopes of Si, only 29 Si has a 
net nuclear spin which may contribute to spectral diffusion by 
flip-flopping with nearby 29 Si. Natural Si contains about 5% 
29 Si or less through isotopic purification. Isotopic purification 
or nuclear polarization will suppress spectral diffusion. 

of our cluster expansion and provide independent verifi- 
cation of its validity; and finally in Sec. IVII we apply 
our theory to obtain numerical results. We conclude in 
Sec. IVIII Five appendixes provide some mathematical 
details and proofs. 

II. BACKGROUND 

In systems of interest to us, a qubit is represented by 
the spin of a localized electron, at low temperature and 
in a strong external magnetic field, with a wave func- 
tion that typically envelopes hundreds of thousands of 
nuclei or more in the surrounding neighborhood of the 
lattice. Two such examples are the shallow donor elec- 
tron state in a phosphorus dopant in silicon- and a local- 
ized electron in a quantum dot in GaAs£ A schematio2& 
of the spectral diffusion process in the former is shown in 
Fig. [J] As discussed in the Introduction, the dominant 
decoherence method is spectral diffusion (SD) caused by 
the dynamic evolution of nuclear spins resulting in effec- 
tive magnetic field fluctuations experienced by the elec- 
tron. The electron spin couples to nuclear spins via hy- 
perfine interaction in the region in which the electron 
wave function is appreciable. The nuclear spins couple 
to each other via dipolar interactions. The large external 
field suppresses processes that do not conserve electron 
or nuclear spin polarization in the direction of the mag- 
netic field. Thus the only relevant nuclear processes are 
flip-flops in which one spin is raised while the other is 
lowered. 

SD is a dephasing decoherence (i.e., a transverse or T 2 - 
type relaxation) process, affecting only the precession of 



the electron spin in the Bloch sphere without changing 
the electron spin component along the magnetic field. It 
thus contributes to the energy-conserving T 2 decoherence 
rather than T% decoherence in which energy is exchanged 
with the bath (Ref. H3 details our definition of T\, T 2 , 
and T 2 *). The T\ time for these systems at low tem- 
perature is known to be much longer than this T 2 time. 
To analyze this decoherence, we consider an ensemble of 
spins, initialized in some direction perpendicular to the 
magnetic field, and observe the decay of the average spin 
over time. In experiments, many sparse donors or dots 
serve as the ensemble, and they are initialized by apply- 
ing a 7r/2-pulse to rotate electron spins polarized along 
the magnetic field axis into the plane perpendicular to 
the magnetic field. For the free induction decay (FID), 
no further pulses are applied and the observed decay is 
largely due to inhomogeneous broadening which is caused 
by the difference in the local magnetic field of each elec- 
tron, causing precession at different rates for different 
electrons. This will provide what is referred to as the T 2 * 
decoherence time. At the very least, this inhomogeneous 
broadening is a result of the random initial distribution of 
nuclear spin states which gives an uncertainty that scales 
as y/N where N is the effective number of nuclei (> 10 5 ) 
influencing the electron significantly. This decoherence 
due to inhomogeneous broadening can be eliminated by 
rcfocusing methods such as the Hahn echo. 24 In such an 
experiment, a 7r-pulse, resonant with electron spins, is ap- 
plied at time t = r in the same direction as the original 
7r/2 pulse and then an echo is observed at t = It. This 
has the effect of reversing any static local field effects and 
allows us to obtain the actual SD decoherence (i.e., the 
T 2 decoherence as opposed to T 2 decoherence) caused 
by effective magnetic field fluctuations. This experiment 
provides a measure of T 2 . The current paper will analyze 
the Hahn echo experiment. More sophisticated pulses, 
such as the Carr-Purcell-Meiboom-Gill sequence^ can 
yield even longer decoherence times and will be analyzed 
in future publications. 

Previous attempts at analyzing this SD decoherence 
have been based on quasiclassical stochastic modeling. 
Herzog and Hahni assigned a phenomenological Gaus- 
sian probability distribution function for the Zeeman 
frequency of the investigated spin without considering 
the dynamics of the nuclear bath. Later, Klauder and 
Anderson* used a Lorentzian distribution function in- 
stead in order to account for a power law time depen- 
dence observed in experiments by Mims and Nassau. 3 
Zhidomirov and Salikho\si devised a more sophisticated 
theory, with a wider range of applicability, in which the 
flip rate of each spin in the bath was characterized by 
Poisson distributions. Very recently de Sousa and Das 
Sarma, 9 in considering spectral diffusion by nuclear spin 
flip-flops, extended this theory to characterize flip-flop 
rates of pairs rather than individual spins within a phe- 
nomenological model. 

In this paper, we present a microscopic theory that 
is based entirely on the quantum mechanics of the sys- 
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tern without resorting to phenomenological distribution 
functions. No Markovian assumption nor any assump- 
tion about the form of the solution was used to obtain 
our results here. Our expansion is entirely based upon a 
density matrix formulation of the problem which assumes 
infinite nuclear spin temperature (T 3> nK) and uses an 
approximate but microscopic Hamiltonian. The problem 
obviously involves too many nuclear spins to solve di- 
rectly using exact Hamiltonian diagonalization; however, 
the cluster expansion method we devise can give succes- 
sive approximations to the exact solution (convergent for 
short times, but often out to the tail of the decay such 
that the full solution is obtained for practical purposes). 
A short report of our results without any theoretical de- 
tails has earlier appeared as a Rapid Communication • 
Our lowest order solution^ was recently reproduced by 
Yao et alm& using an entirely different approach, thus 
providing an independent validation of our theoretical 
approach. 



III. FORMULATION OF THE PROBLEM 

A. Free evolution Hamiltonian 

The free evolution Hamiltonian for the spectral diffu- 
sion problem, considering one localized electron spin, S, 
and a nuclear spin bath, J„, in an external magnetic field 
B (taken to be along the z direction), is given by 



n = n 



Ze 



H 



Zn 



Tt + TL 



(1) 



where the first two terms are due to the Zeeman ener- 
gies of the electron and nuclei respectively, TL gives the 
electron-nuclei coupling, and TL B contains the internu- 
clear coupling. 

The Zeeman energy contributions, arising from the ex- 
ternal magnetic field, are given by 



tl z 



= lsBS z , 



(2) 
(3) 



where 75 and "f n are the gyromagnetic ratios of the elec- 
tron spin and nuclear spins, respectively, and B is the 
external magnetic field defined to point in the z direc- 
tion. 

The electron-nuclei coupling is given by 



H A = ^,A n (I n .S). 



(4) 



A n gives the magnetic coupling between the electron spin 
and each nuclear spin. It is dominated by the hypcrfinc 
coupling for nuclei that contribute significantly to SD. 
Since 75 = 1.76 x 10 7 (s G) -1 is typically four orders of 
magnitude larger than 7„, flip-flop processes between the 
electron and a nuclear spin are strongly suppressed at low 
temperatures (i.e., no phonons) by energy conservation 



in a strong magnetic field^ This direct hyperfine inter- 
action leads to quantitatively important effects at zero 
magnetic field^S but at the moderate or high magnetic 
fields required for spin resonance measurements it only 
contributes a small visibility decay^ We will therefore 
disregard the nonsecular part of the electron-nuclei cou- 
pling, leaving us with 



TL A ~ ^ ' A n I nz S z 



(5) 



which we will use as the definition for TL A throughout the 
rest of this paper. 

The internuclear coupling due to the dipolar interac- 
tion is given by^i 



TL 



<T_/D 



B 

nm 1 



R 3 

' n in 



R 5 



(6) 



,(7) 



where R nm is the vector joining nuclei n and m. This 
can be expanded into a form containing only operators of 
the type I + , i_ , or I z M- We will again invoke the energy 
conservation argument which, because of the strong ex- 
ternal magnetic field, allows us to neglect any term that 
changes the total Zeeman energy of the nuclei. This will 
leave us with the following secular contribution: 



Ti„ 



u n in n ~ lm)In+Im- ~ ^l-nzlmz), (8) 



1 1 - 3 COS 2 I 



(9) 



where 9 nm is the angle of R nm relative to the mag- 
netic field direction. Equation JSJ will be used for H^m 
throughout the rest of this paper. Note that the flip-flop 
interaction between nuclei with different gyromagnetic 
ratios is suppressed by Zeeman energy conservation in the 
same way that the nonsecular part of the dipolar inter- 
action is suppressed. This occurs, for example, in GaAs 
because the two isotopes of Ga and the one isotope of As 
that are present have significantly different gyromagnetic 
ratios. This approximation is applicable if any such dif- 
ferences in 7„ are on the order of |7„| itself and thus such 
energy changes would be large compared to other ener- 
gies of the problem in the large external magnetic field 
being considered in this work. Giving some typical num- 
bers, max(|A n |) ~ 10 6 s _1 , max(|6„ m |) ~ 10 2 s _1 , and 
for B = 10 T, \ ln B\ - 10 8 s" 1 and \ ls B\ - 10 12 s" 1 . 
The Zeeman energies are thus much larger than the other 
energies in the problem. In fact, our model of neglecting 
terms that change the total Zeeman energy of the nuclei 
turns out to be reasonably valid for B > 0.1 T or so 
for Si:P or GaAs systems of interest to us. SD is there- 
fore the dominant spin decoherence channel for localized 
electrons in semiconductors except at zero or very small 
(< 0.1 T) magnetic fields. 
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B. Hahn echo 

The Hahn echo experiment consists of preparing the 
electron spin in a state perpendicular to the magnetic 
field, allowing free evolution for a time r, applying an 
electron spin resonant 7r-pulse about an axis perpendicu- 
lar to the magnetic field (for our theoretical purpose, the 
direction of this axis is not important beyond being per- 
pendicular to the B-field), and then allowing free evolu- 
tion again for time r. The echo envelope is the magnitude 
of the resulting ensemble spin average at time t = It. 

Mathematically, the echo envelope is the magnitude 
of the expectation value of the spin multiplied by 2 for 
normalization (since the electron has a spin of 1/2). Be- 
cause we are only concerned with the component of the 
spin perpindicular to the magnetic field, it is also the 
magnitude of the normalized complex in-plane magneti- 
zation: 

v e {t) = 2((S x )+i{S y )) (10) 
= 2 Tr {(S x +iS y )p(r)}. (11) 

The density matrix of the system, p(r), is given by 

p(r) = U(t) Pq U\t), (12) 

with the evolution operator 

U(t) = e- mT a x ^ m \ (13) 

where TL is the free evolution Hamiltonian given by 
Eq. QJ. We have arbitrarily chosen the pulse axis to 
be in the x direction. This pulse is represented by the 
Pauli matrix, <j xe — 2S X with the e subscript denoting 
that it is an electron spin operator. 

Our initial density matrix, po, will be given by a prod- 
uct state of the electron spin, in a direction perpendicular 
to the magnetic field, and the nuclear spins in a thermal 
state: 

Po = \xl){xl\® e ~ U ^ BT , (14) 
lx"> = ^(|0 e ) + e^|le», (15) 

where Ti n = H Zn + TL B and M is its partition function 
(M rj 2^ for T 3> nK, where TV is the number of nuclear 
spins). The initial angle of the spin, <f>, will contribute 
an irrelevant phase factor to the complex in-plane mag- 
netization. Dropping this irrelevant phase (since we are 
interested in the magnitude), after a few manipulations, 
shown in Appendix 1X1 Eq. (|llfl becomes 

Mr) = ^ Tr [U-U + e- Hn / ksT ulul}, (16) 

where 

U± = e- m±T (17) 



are evolution operators under the effective Hamiltonians 

n^m n 

which describe nuclear evolution under the effect of an 
electron spin up (H + ) or down (Ti~) with the exter- 
nal magnetic field dependence now canceled out. Al- 
though the U± evolution operators are independent of 
B, it should be noted that a strong external magnetic 
field is required to justify the approximations of Eqs. (jSJ) 
and © . The trace in Eq. (|16fl is taken over nuclear spin 
states only. 

For our purposes we only consider high nuclear tem- 
peratures (T 3> nK) such that we can make the following 
approximation: 

v E (r) = ^Tr{U-U + Ulut}. (19) 

The high nuclear temperature (and the low electron and 
phonon temperature) approximation made throughout 
this work is well-valid in the solid-state quantum com- 
puter architectures of our interest, where typical experi- 
ments would be carried out around T ~ 100 mK which 
is an extremely high temperature for the nuclear bath 
dynamics and an extremely low temperature for phonon 
excitations. 



IV. CLUSTER METHOD 

Our cluster method provides a strategy to compute 
the Hahn echo, Ve(t) [Eq. JTJJ], in successive orders of 
accuracy which may or may not converge for a specific 
application (naturally the systems reported in this pa- 
per do converge, for practicle purposes, as we will show). 
We refer to these successive orders, loosely, as an "expan- 
sion" although we do not wish to suggest that corrections 
are additive (as in familiar expansions such as a Taylor 
series). In Sec. IIV Al we provide a conceptual descrip- 
tion of this cluster expansion. Section IIV Bl supplies the 
mathematical formalism and practical implementation of 
this expansion. 

A. Conceptual cluster expansion 

Consider independent, simultaneous nuclear "pro- 
cesses" that may contribute to the decay of the Hahn 
echo. For example, a process may involve a pair of nuclei 
flip-flopping [Fig. |2fa)] which results in fluctuations of 
the effective magnetic field seen by the electron spin, or 
it may involve three nuclei interdependently [Fig. 12Tb)]. 
etc. The dynamics of such a process results from the 
local coupling between nuclei (with coupling constants 
{b nm }), and hyperfine coupling to the electron (with cou- 
pling constants {A„}). Any number of these processes 
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FIG. 2: Some possible nuclear "processes" where the num- 
bered two-sided arrows represent a sequence of flip-flops be- 
tween pairs of nuclei, (a) depicts a two-nuclei process and (b) 
depicts a three-nuclei process. 



may occur "simultaneously" as long as they involve dis- 
joint sets of nuclei and are thus independent of each other 
(processes that share a nucleus are not independent and 
would have to be combined into a larger process). 

Using this (not yet well-defined) concept of nuclear 
processes, the cluster expansion may be described, ide- 
ally, as follows. The cluster expansion will include pro- 
cesses that involve a successively increasing number of 
nuclei. Because an isolated nucleus in our model does 
not contribute to spectral diffusion, at the lowest nontriv- 
ial order we include any simultaneous processes involving 
two nuclei (pairs). Thus at this lowest order we can in- 
volve any number of pair processes together as long as 
the pairs do not overlap (i.e., involve the same nucleus). 
At the next order, we will additionally include processes 
that involve three nuclei. Next, we include four nuclei 
processes which cannot be decomposed into two pair pro- 
cesses (these would already have been included) . To sum- 
marize, let us say that the k th order of the expansion will 
include processes of up to k nuclei. 

Because all processes involving a given number of nu- 
clei are included at each order of this expansion, and be- 
cause these processes are independent (proven formally 
in Sec. IIV Bl and the related appendix), rather than work- 
ing with individual processes, we can work with contri- 
butions due to each given "cluster" of nuclei (for now, 
simply defined as a set of nuclei); such a "cluster contri- 
bution" includes contributions from all of the processes 
involving all nuclei in that cluster in an interdependent 
way (i.e., not separable into independent sub-processes). 
Thus we may say that the k th order of the expansion 
includes contributions from clusters up to size k. These 
"contributions" are not necessarily additive in the solu- 
tion because we must account for simultanous but inde- 
pendent processes (from disjoint clusters). The idea is 
simply to include the possibility of interdependent pro- 
cesses involving clusters of successively increasing size. 

We deliberately use the word "cluster" to imply prox- 
imity between the members of the set of nuclei involved 
in interdependent processes. In fact, a near neighbor 
approximation, in which the constituent nuclei of a con- 
tributing cluster must be in the same neighborhood, is 
justified by the 1/R^ m dependence of the internuclear 
coupling constant [Eq. @]. Consider a near (not nec- 
essarily nearest) neighbor approximation with an ad- 




FIG. 3: L is loosely defined as the average number of neigh- 
bors in a near neighbor approximation that converges to the 
exact answer. For example, we can include neighbors up 
to a distance ro away such that increasing this maximum 
neighbor distance in the near neighbor approximation (where 
non-neighbor interactions are neglected) does not significantly 
change the solution. 



justable parameter r such that we ignore interactions be- 
tween nuclei that are a further apart than r. If a near 
neighbor approximation is applicable, the Hahn echo so- 
lution in this approximation (in principle, whether or not 
it is feasible to compute) will converge with an acceptable 
level of accuracy at some finite value of r much smaller 
than the system size. Let us define ro to be the value of 
r in which this acceptable convergence is achieved. Let 
L be the number of nuclei within a range of ro from any 
nucleus, on average, as shown in Fig. |3| Applying this 
near neighbor approximation to our cluster expansion, L 
determines the way in which the number of contribut- 
ing clusters scales with cluster size. This has important 
implications for the convergence of the cluster expansion. 

To be specific, the convergence of the cluster expansion 
depends upon two factors. The first is how the number 
of contributing clusters scales with cluster size (which re- 
lates to L as we have already said). The second is how 
the average contribution of clusters scales with cluster 
size. Clusters only contribute via interdependent pro- 
cesses; thus the set of nuclei in a contributing cluster 
must form a connected graph where edges in the graph 
connect neighbors [Fig. Ufa)]. When counting the num- 
ber of clusters of a given size, we have N sites to choose 
from for the first nucleus, but there are only O(L) pos- 
sibilities (roughly) for each additional nucleus because it 
must neighbor one of the previous choices. This sim- 
ple analysis does not compensate for overcounting due to 
permuting labels and other such details, but it provides 
the correct scaling in terms of N and L; that is, there 
are 0(NL k ~ 1 ) contributing clusters of size k. Our first 
scaling factor is then L since the number of clusters as a 
function of cluster size scales in powers of L. The other 
scaling factor will rely upon some perturbation theory to 
describe how cluster contributions themselves scale with 
an increase in cluster size. Section will discuss two 
perturbation theories that may be applicable [for short 
times, r <C max (t nm ) ], but for now we will use A to 
represent some perturbation parameter and assume that 
a cluster contribution of size k will scale as 0{X k ). Thus A 
is assigned as the other scaling factor, and we may loosely 
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a) valid invalid 




b) b 1+ b, ; b +s b 5S (b i7 ) 2 (b 36 ) 2 




FIG. 4: (Color online) (a) The set of nuclei in a contribut- 
ing cluster must form a connected graph. Edges represent 
neighbor connections. The set of blue nuclei on the left form 
a valid cluster. The set of green nuclei on the upper right 
are not fully connected so they do not form a valid cluster, 
(b) A set of bnm factors in a term of an expansion of Ve(t) 
determines a set of disjoint clusters (the connected subgraphs 
formed from b„ m edges). 



argue that we expect the cluster expansion to converge 
when XL <C 1 because the total contribution from clus- 
ters then decreases as we increase in cluster size. This 
reasoning will become more rigorous when we go on to 
explain how we implement this cluster expansion. 



B. Mathematically formalized cluster expansion 



In Sec. IIV Al we gave a rough, conceptual description 
of our cluster expansion to guide the reader's intuition 
and present some basic ideas. At this point, we will de- 
velop the rigorous mathematical formalism that relates 
the idea of simultaneous, independent nuclear processes 
contributing to the Hahn echo directly to our mathemat- 
ical expression [Eq. Q19[l] for the Hahn echo. We will 
decompose Eq. (I19fl into a sum of products of cluster 
contributions. Each cluster contribution will effectively 
contain the sum of contributions from all processes in- 
volving, inseparably (i.e., interdependently) , all nuclei in 
the cluster. Such a decomposition requires that processes 
involving disjoint sets of nuclei are truly independent and 
interchangeable. This requirement is met by proving, as 
we shall, that a cluster contribution is independent of 
anything outside of the cluster. 

These cluster contributions need not be computed by 
analyzing the possible "processes" involving each set of 
nuclei. Instead, the decomposition of Eq. l(T9")l into clus- 
ter contributions will be used recursively to define these 
cluster contributions; this is shown in Sec. IIVB II With 
these cluster contributions concretely defined, we then 
discuss, in Sec. IIV B 21 how we mathematically define 
the ideal cluster expansion that we have conceptually de- 
scribed. This ideal expansion is useful for understanding 
some basic ideas, but in order to practically perform cal- 
culations on large systems, some further approximation 
techniques must be used. This practical implementation 
of the cluster expansion is explained in Sec. IIV B 31 



1. Decomposing into cluster contributions 

Consider expanding ve(t) [Eq. (|T^|) ] into a sum of 
products with respect to internuclear coupling, having all 
bnm constants factorable from each term. For example, 
such an expansion could be made by Taylor expanding 
the exponentials of U± [Eq. l|17fl] and then distributing 
through these sums. Each term in such an infinite ex- 
pansion has a set of b nm factors which determine a set 
of involved clusters. In the language of graph theory, the 
bnm factors may be represented by edges (between nodes 
n and to); then the clusters are the sets of nuclei in each 
connected subgraph (each of these being involved in an 
interdependent process). Figure Efb) illustrates an ex- 
ample of this. In this way we begin to relate our concept 
of clusters of interdependent processes to Eq. . 

Each term in such a sum of products expansion will 
involve some disjoint set of clusters. Let us assume that 
processes involving disjoint clusters are truly indepen- 
dent. For a particular cluster, C, we should then be able 
to factor out a unique cluster contribution from those 
terms which involve cluster C . This is the basis for a de- 
composition of ve{t) into cluster contributions. It will 
be useful to generalize this by defining vs{t) to be the 
solution to the Hahn echo, ve(t) [Eq. (|19|l ]. when only 
including the nuclei in some set S. We will use v' c (t) to 
represent the contribution from cluster C. Appendix iBl 
proves that the cluster contributions are independent of 
each other in a way that allows for the following decom- 
position of vs{t): 



Vs(t) 



E 

{Ci} disjoint, 
d ^ 0, d C S 



i+ E IRw. 



(20) 



(21) 



{Ci}#0 disjoint, i 
C, jt 0, C, C S 



where the summation of Eq. I|2(J|) is over all possible sets, 
{Ci}, of disjoint nonempty clusters, Ci, each of which is 
contained in or equal to S. In other words, it iterates 
over all possible ways of dividing any part of S into dis- 
joint clusters as depicted in Fig. The product is over 
all clusters in each set. Despite the index, i, the or- 
der is irrelevent and permutations do not count as dis- 
tinct cases. Extracting the trivial {Ci} = term yields 
Eq. i|21|) . shown explicitly to avoid confusion or ambigu- 
ity. The unique existence of such a decomposition follows 
from the fact that any v' c {r) must be well-defined inde- 
pendent of any nuclei outside of C which is proven in 
Appendix IBl 

We can use Eq. (|20|) itself to obtain an unambiguous 
expression for any v' c (r). We do this by applying Eq. lt2U|) 
to the case in which S — C and pulling out the term, from 
the summation, in which {Ci} = {C} leaving only sets in 
which all C % 7^ C: 



v c(t) = v' c (t) 



E IRw 



(22) 



{Ci} disjoint, 
Ci^0, d C C 
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FIG. 5: Set of all possible sets, {Ci}, of disjoint contributing 
clusters contained in a set, 5, of four nuclei as an example. 
Contributing clusters are of size 2 or greater (a single nucleus 
gives no contribution on its own). The cases on the left in- 
volve 2-nuclei, middle ones involve 3-nuclei, and the ones on 
the right are the trivial cases of {Ci} — or {Ci} — {S}. Such 
possibilities are iterated over in the summation of Eq. 1201 . 
Equation 121H removes the {Ci} = case and Eq. 1231 re- 
moves the {Ci} = {5} = {C} case from their respective sum- 
mations. 



so that 



v c( T ) = v c{t) 



E 

{Ci} disjoint, 
Ci^0, C x C c 



(23) 



Equation ill' Ml) provides a recursive definition of a cluster 
contribution. Starting with the computation of vc(t), 
which is feasible to calculate directly for small clusters, 
one must subtract terms that involve multiple indepen- 
dent processes and processes that do not involve all of 
the nuclei in C. It is a direct consequence of the decom- 
position given by Eq. 1(20) 1 . 

To ensure that Eq. ill' /it is well-understood, we show 
more explicit results for clusters of size one through four. 
Because a single isolated nucleus cannot contribute to 
spectral diffusion, v' c (r) = vc t (r) — 1 = when \C\\ = 1. 
It follows that for 2-clusters, v' C2 (r) = vc 2 (r) — 1 (with 
\C2\ = 2), having no contributing proper subclusters. For 
3-clusters, we must subtract off contributions from con- 
tained pairs: 



(r) = v C3 (t) - 1 - 



E 

C 2 C C 3 , 
|C 2 |=2 



(24) 



For 4-clusters, we must also subtract off contributions 
from contained 3-clusters and the products of contribu- 
tions from contained disjoint pairs: 

v' Ci (r) = v Ci (r)-l- £ v' C2 (r)- £ v'^r) 



C 2 C Ci, 
|C 2 |=2 



c 3 c c 4 , 

|C 3 |=3 



E ^W«bW« 



(25) 



A U B = C 4 , 
\A\=\B\=2 



The factor of one-half in the last term is needed to com- 
pensate for the fact that A and B may be swapped in the 



summation; it is only a consequence of the notation used 
here (where A and B are interchangeable labels). 



2. Ideal cluster expansion 

We are now able to compute cluster contributions to be 
used in the evaluation of our cluster expansion. Revising 
Eq. (|20() slightly, we may write the following expression 



for the ideal cluster expansion up to k order 
4\r)= £ l[v' Ci (r). 

{Ci} disjoint, i 
C ( #0, |C,|<fe 



(26) 



In order to estimate the error of the k th order of the ex- 
pansion, we can compare it with the (fc + l) th order which 
must include contributions from fc+1 sized clusters. One 
way to convert (r) into (t) is to add additional 

terms to the sum in which we replace any fc-cluster con- 
tribution of an existing term with any (k+ l)-cluster con- 
tribution generated by adding one neighboring nucleus to 
the original fc-clustcr. In doing so, a replacement must be 
made because the original fc-cluster becomes disqualified 
when we introduce the new (k+ l)-cluster which contains 
it (due to the requirement that the clusters be disjoint). 
This approach will account for all new sets of {Ci} con- 
taining a (k + l)-cluster (since any (k + l)-cluster can be 
made by adding a nucleus to a fc-cluster); however, cases 
will be overcounted because many fc-clusters can be used 
to build the same (fc -I- l)-clustcr. This is unimportant 
because our goal now is to estimate the error of (r) 
relative to v^ +1 \t) and overestimating this error is just 
as good. Proceeding along these lines, we first separate 
out the fc-cluster contributions: 



,(*0 



(r) 



E 

{Ci.T>j} disjoint, 
0<\Ci\<k, |X>,-|=fc 



IKwIKa)- ( 27 ) 



This performs the same summation over sets of disjoint 
clusters as in Eq. (|26() except that we label fc-clusters as 
T>j and the smaller clusters as Ci. With these fc-clusters 

now set apart, we can estimate the error of v^\t) rela- 
tive to v^ +1 \t) by noting that the sum of all (fc + 1)- 
cluster contribution replacements of v' v . (r) are roughly 
O(XL) x Recall that A was introduced as a 

perturbation parameter such that a cluster contribution 
of size fc scales as 0(X k ), and L is the average number 
of neighbors so that there are, roughly speaking, O(L) 
(fc + l)-clusters that may be built out of one fc-cluster. 
Thus 



v 



(fe+i) 

E 



(r) = 



E IRMx 



{Ci,Vj} disjoint, i 
0<\d\<k, \Vj\=k 

Y[v' v .{t) [1 + 0(XL)}. 



(28) 
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FIG. 6: One possible combination of simultaneously included 
pair contributions. The red/dark circles are the nuclei whose 
processes are being considered. 



If we explicity include these (k + l)-clusters, they would 
have relative corrections of O(XL) to account for (k + 2)- 
clusters and so forth. This provides a more rigorous argu- 
ment for our previous assertion that the cluster expansion 
converges when XL <C 1. 



3. Practical implementation of the cluster expansion 

Equation i|26|) directly implements the conceptual clus- 
ter expansion as described in Sec. IIV A*l (the inclusion of 
contributions from all clusters up to size k); however, it 
is impractical for calculating results in large systems. At 
the lowest nontrivial order, we would need to sum over 
all possible products of disjoint pair contributions; for 
example, Fig. [H] depicts one such combination of disjoint 
pairs. The number of such possibilities likely grows ex- 
ponentially with the problem size (we have not bothered 
to prove this rigorously or look it up, but certainly the 
scaling with problem size is horrendous). However, we 
can effectively obtain all possible combinations by mak- 
ing products of the form FJ C [1 + v' c {t)\. Distributing 
through a given factor yields the possibility of excluding, 
via the 1 term, or including, via the v' c (r) term, that clus- 
ter. Therefore such a product gives the sum of all pos- 
sible combinations of simultaneous cluster processes (for 
the clusters included in the product). Unfortunately, this 
will yield combinations that involve overlapping clusters 
(that are therefore not independent). These overlapping 
clusters will introduce an error that, in principle, may be 
corrected in successive orders of an approximation. 

With this approach, the lowest nontrivial order of the 
expansion may be implemented with 



*4 2) (t)« n 

|C 2 |=2 



(29) 



producing all combinations of pair contributions along 
with some extraneous terms, such as overlapping pairs 
as depicted in Fig.[7fb). For a moment, let us disregard 
these erroneous terms and consider the consequence of 
this approximation. If we take the logarithm of both 
sides, we can convert the product on the right-hand side 



of Eq. I|29[) into a convenient sum: 

« ln(l + <4(T)) (30) 

|C 2 |=2 

« E ^(r)[l + 0{v' C2 (r))], (31) 



\C 2 \=2 



where Eq. I|31|) follows from the Taylor expansion of 
ln(l + w^ 2 (r)J for |vc 2 ( T )| ^ 1 which we will shortly 
justify in a self-consistent way. If we assume that \v' c (r)| 
is small for all (or most) of the C 2 pairs, then 



4V) 

£ 2 (r) 



« exp[£ 2 (r)], 

= E <4M- 



(32) 
(33) 



|C 2 |=2 



It is easy to show, from Eq. I|19l) . that vs(t) is al- 
ways real and that —1 < wc 2 ( T ) < 1> therefore, —2 < 
[v' c% {t) =v C2 {t) - 1] < 0. Because «^(t = 0) = 
vc 2 (t — 0) — 1 = 0, £2(7" = 0) = and becomes in- 
creasingly negative (initially at the very least) as r is in- 
creased. For a large system, we expect S 2 (r) to decrease 



monotonically to a value that is 0{— N) [i.e. 



,(t) 



have become random] so that Eq. I|32|l exhibits a de- 
cay form. The interesting part of the decay occurs when 
^2(7") > —1 so that v^\t) > er 1 . When this is true, 
the average pair contribution will be, at most, 0(1/N), 
self-consistently justifying the approximation of Eq. I|32f) 
relative to Eq. I|31|l when N is large (as it is for systems of 
interest). Increasing r much beyond this point will bring 



us to the tail of the decay in which ve (t) 



,(2) 



(r) « 1. 



To state this in a physically intuitive way, the decoher- 
ence of spectral diffusion is caused by many nuclei collec- 
tively such that each potentially flip- flopping nuclear pair 
contributes only a small amount to the overall dephasing 
before coherence is completely lost. 

For practical purposes [when u_e(t) is not terribly 
small], we thus regard each cluster contribution to be 
(D(\/N). Now let us discuss the extraneous overlapping 
pairs of Eq. l(2l?|) that we have thus far disregarded. We 
can now think of these cases, and their corrections, in 
orders of 1 /N which increase with the number of overlap- 
ping clusters. The lowest order correction will therefore 
remove cases of two pairs that overlap with each other. 
For any given pair, there are O(L) pairs that can over- 
lap with it, each of which has a contribution of 0(1/N) 
as discussed above. Therefore in estimating the error of 
Eq. I|29|l we may write 



*4 2) w= n (i+vkw 



ic 2 =2 



1 



(34) 



Carrying this error over to Eq. (|32|l and including its 
inherent 0(1/ N) error, we have 



4 2) (r) =exp( E 2 (r) 



l + O 



N 



(35) 
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FIG. 7: The practical implementation of the cluster expan- 
sion approximates the ideal cluster expansion up to errors 
resulting from overlapping clusters. At the lowest order, such 
errors involve overlapping pairs, (a) A single pair multiplied 
by itself (i.e., a pair overlapping itself) and (b) two pairs over- 
lapping by sharing a nucleus. 



Because a cluster contribution scales in orders of A 
as we increase the cluster size, this approach may be 
used for higher order cluster contributions provided that 
A <C N (typically, A <C 1 where the cluster expansion is 
applicable). Taking either A <C 1 or A ~ 1, we may write, 
as an extension of the above approach to higher orders, 



Sj(r) = £ v' c (r). 

\C\=3 



1 + 



l + L 

N 



(36) 
(37) 



Note that £fc(r) ~ T,k-i(r) x O(XL), since there are 
roughly O(L) times as many fc-clusters as (k — l)-clusters 
and on average each fc-cluster contribution, by the def- 
inition of A, is 0(A) times that of the average (k — 1)- 
cluster. With this in mind, we see that, under the cluster 
expansion, 1il(ve(t)) is effectively expanded, additively, 
in powers of (XL). 

In addition to the expansion in cluster size, we may also 
successively correct for the 0(l/N) errors of overlapping 
clusters. This is done by starting with the smallest num- 
ber of overlapping clusters of the smallest sizes; that is, 
start with the case of two overlapping pairs (Fig. 01 . Each 
additional cluster included in the set of overlapping clus- 
ters being considered will multiply 0(1/N) to the cor- 
rection, and each additional nucleus added to any cluster 
will multiply 0(AL) to the correction. For our purposes, 
we will only consider the correction for two overlapping 
pairs as a check to verify that the approximation made 
in Eq. (J2HJ) is valid. 

There are two cases to consider for this lowest order 
correction of overlapping clusters: the same pair multi- 
plied by itself [Fig. 0a)] which was introduced by the 
approximation of Eq. I|31() . and two different pairs that 
overlap [Fig. 0b)] which originates from Eq. J2HI- These 
cases are, respectively, elimated, to lowest order (two 
pairs and only two pairs that overlap), by adding the 



following to Eq. ffij} : 

S 2(r) = ~\ E KM]' 



so that 

\n(yf{r)) 



\C 2 \=2 

\ E Va(tW b (t), 



(38) 
(39) 



\A\JB\=3, 
|X|=|B|=2 



£ 2 (t) + £*(t) + £*(t) 
k 



(40) 



E s .( 7 

J=3 



1 + 



1 



N 



Exponentiating Eq. I|40|) then expanding and distributing 
this exponential into a sum of products form will yield 
the sum of all products of disjoint cluster contributions, 
as in Eq. I|2t)l) , plus extraneous terms of overlapping clus- 
ters. However, all cases of only two pairs overlapping 
with each other (including a pair multiplied by itself) will 
be removed as a result adding in S^t) and S|(r). There 
will remain higher order errors with more than two over- 
lapping clusters or overlapping clusters larger than pairs; 
in fact, additional higher order errors are introduced by 
the £j (r) and £3(7") corrections itself. For this reason, 
it is difficult to derive higher order corrections (you must 
correct errors introduced by lower order corrections). We 
can, however, regard this lowest order correction as an 
estimate of the error caused by these extraneous overlap- 
ping clusters: 

k 

la(4 fc) to) = J2^(t)+0(Z*(t)), (41) 
£^(r) + £*(r). 



£*(t) 



(42) 



Note that £^(r) and £3(7") are both < and there- 
fore add constructively (otherwise we would want to take 
absolute values in order to estimate the error conserva- 
tively). Fortunately, calculations of £*(r) indicate that it 
is a minor correction for practical purposes. Such calcu- 
lations verify the argument that these are 0(1/ N) errors 
[at least for practical values of r for which Ve(t) > e^ 1 ]- 



C. Cluster expansion in summary 

The cluster expansion method that we have developed 
in this section is very powerful and very general. The 
disjoint cluster decomposition [Eq. I|20(l ] could be used 
to take the trace of any evolution operators described by 
Hamiltonians with pairwise (or even higher order) inter- 
actions. This decomposition may then be used to form 
an expansion [Eq. I|26|l ] that converges when the sum 
of cluster contributions decreases with cluster size (i.e., 
XL < 1). In order to practically compute this expan- 
sion for a large system, we need to use approximations 
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such as Eq. (|34|) or Eq. (|36|l which have the additional 
requirement that each cluster contribution be small, e.g., 
0(1/N), so that extraneous overlapping clusters arising 
from these approximations are small. This is, in princi- 
ple, a formally exact, systematic expansion and its con- 
vergence may be tested by comparing Sj(r) for at least 
j = 2, 3, and 4 as well as S*(r). It is important to com- 
pute £4(7") as well as S 3 (r) because a symmetry of the 
problem eliminates all odd orders of A as we will show for 
both perturbation theories in Sec. [VJ therefore, 3 cluster- 
contributions are actually 0(A 4 ). 

We conclude this section by remarking that, besides 
being elegant and useful for understanding the expan- 
sion, the natural logarithm form of the Hahn echo given 
by Eq. (|41H has the advantage that it is convenient to 
compute Sj(r) and S*(r) using Monte Carlo techniques. 
Rather than computing the full sum, randomly selected 
terms may be sampled and averaged in order to obtain 
an estimate for each sum. This can save a lot of com- 
putation time and makes this method powerful for large, 
complicated systems. 

V. DUAL PERTURBATION THEORIES 

The convergence of the cluster expansion, described in 
the previous section, depends upon how cluster contribu- 
tions scale with cluster size. In that section, we surmized 
that a perturbation theory could be invoked such that 
cluster contributions scale with cluster size in orders of 
some perturbation parameter, A; that is, the contribution 
of a cluster of size k is 0(X k ). In this way, the cluster 
expansion relies upon an underlying perturbation theory 
to justify or explain its convergence. Alternatively, one 
may regard the cluster expansion as a means to extend 
convergence of a perturbative expansion to large systems 
(i.e., even when N 3> A). In any case, in order to form a 
connection between a perturbation theory and the cluster 
expansion, the contribution of a cluster of size k must be 
shown to have a minimum perturbative order of X k (i.e., 
lower orders cancel out). In this section, we present two 
such perturbation theories which are applicable in com- 
plimentary regimes (and thus, "dual"). 

A. Two theories in complementary regimes 

Although otherwise complementary both perturbation 
theories require that r <C max(& nm ) . However, in 
problems we have considered td <C max(6 nm ) where 
td is the decay time, and therefore this constraint has 
no practical consequence. The first perturbation theory 
we will present is an expansion in orders of r. For this 
theory, we can effectively assign A = max(6„ m )r as the 
perturbation parameter used in Sec. IIVI It is appar- 
ent that r <C max(6„ m ) _1 for A to be small. Besides 
this requirement, we will see that this perturbative ex- 
pansion is generally convergent in the regime in which 



\A n — A m \ <C b nm . The second perturbation theory 
which we call the dipolar perturbation technique, will 
treat H B of Eq. JSJl as a perturbation to the free evolu- 
tion Hamiltonian. For this theory, A ~ b nm /\A n — A rn \ 
and is therefore convergent in the opposite regime as the 
r-expansion with respect to \A n — A m \ versus b nm ; that 
is, b nm \A n — A m \. 

For convenience, we define 



so that the r-expansion perturbation theory is said to 
be applicable in the c nm <C 1 regime while the dipolar 
perturbation is applicable in the c nm 3> 1 regime. The 
factor of 4 in Eq. H43|l really only makes sense in the 
context of spin-1/2 nuclei, but it is irrelevant for the cur- 
rent discussion. Nuclei near the donor or center of the 
dot, where the electron wave function is strongest so that 
hyperfine coupling dominates the internuclear coupling 
[max(^4 n )/ max(6 nm ) ~ 10 3 ], can typically be classified 
in the c nm 3> 1 regime. Nuclei far from the donor or 
center of the dot can typically be classified in the oppo- 
site regime. A consequence of our perturbative theories 
is that the extreme cases will give negligible contribu- 
tions (i.e., there is a "frozen" core where nuclear flip-flops 
are suppressed by the strong hyperfine coupling, and the 
electron's interactions with distant nuclei are too weak to 
be of consequence). In fact, it is arguable that clusters 
of the c nm ~ 1 regime give the strongest contributions 
which seems to contradict our supposition that cither 
of these perturbation theories are applicable. However, 
time plays a role that shifts the balance to the c nm > 1 
side enough to make the dipolar perturbation theory par- 
ticularly applicable in a way that causes the cluster ex- 
pansion to converge. This is because the c nm > 1 regime 
implies larger energy scales (with strong hyperfine cou- 
pling than the converse) which causes such cluster contri- 
butions to operate on shorter time scales and dominate 
the Hahn echo decay. We would not be satisfied with 
the face value of this argument, but we have performed a 
number of numerical calculations that support this con- 
clusion. In particular, we have compared the cluster ex- 
pansion performed in two ways: one in which cluster con- 
tributions are computed exactly, and the other in which 
these contributions are approximated by dipolar pertur- 
bation theory in the lowest order. In order to give con- 
vergent results for the latter calculation, a lower bound 
Cnm cutoff must be imposed. However, for a broad range 
of c nm cutoff values, the two calculations agree very well. 
This implies that dipolar perturbation theory completely 
accounts for the Hahn echo decay. The main importance 
of the r-expansion perturbation theory is to justify the 
c-nm cutoff by the argument that c nrn -C 1 clusters are 
negligible (although it will have other purposes as we 
discuss shortly). 

We emphasize that the convergence of the cluster ex- 
pansion is ultimately proven by performing the calcula- 
tions and comparing the different orders of the expansion 
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(as we do). The dual perturbation theories serve to ex- 
plain this convergence. 



B. Expansion in tau 

Before we devised the cluster method of solving the 
spectral diffusion problem, we attempted to obtain a so- 
lution, at short times, by performing an expansion in r 
directly. This met with failure by not converging (for the 
Si:P system which we tested) even at very short times. 
We now understand that the reason for the failure of the 
direct r-expansion is a combination of the large system 
size and the difference in the energy (time) scales of the 
hyperfine coupling and dipolar coupling. In the context 
of the cluster expansion, however, a perturbation in or- 
ders of r can be applicable. The cluster method deals 
with large systems correctly and can allow underlying 
perturbation theories to exhibit themselves without be- 
ing drowned out by a large N 3> A -1 . 

In Sec. IV B II we will describe the techniques that 
we developed in order to analyze and compute this r- 
expansion. These techniques are shown (in Appendices) 
as they may be adapted to other problems for which they 
are more suitable. In Sec. IVB2I we show that, in the 
regime in which the expansion is convergent for small 
clusters, this provides an underlying perturbation the- 
ory to justify the convergence of the cluster expansion 
(i.e., providing the A used in Sec. IIV|) . Because of the 
difference in the energy (time) scales of the hyperfine 
coupling and dipolar coupling, this perturbation theory 
is not as relevant as the dipolar perturbation described 
in Sec. IV Cl however, we will see that consequences of the 
r-expansion are exhibited in GaAs systems fSec. IVIB)l . 
Also, Fig. 1131 of section I VII provides independent verifi- 
cation of the cluster expansion by using this expansion 
directly (not in the context of the cluster expansion) for 
an artificial system in which it is applicable. 

1. Algebraic expansion 

An expansion in r is performed by expanding the ex- 
ponentials of U±(t) [Eq. fTTjl ] in Eq. (fT§|l and collecting 
all distributed terms of the desired order in r. One can 
additionally expand e~ u l kBT , from Eq. Q16)). in powers 
of T" 1 but in the remaining discussion we will be taking 
the T — > oo limit. We developed a symbolic manipu- 
lation computer program that gives an exact algebraic 
expression for any expansion order (although computa- 
tion time limits the expansions that are feasible). The 
program uses its own algebraic manipulation library cre- 
ated precisely for this purpose. AppendixOdescribes the 
procedure used by the program to obtain this expression. 

The computer program's results for 1=1/2 were com- 
pared to direct hand calculations (as a check) for up to 
sixth order in r. Actually, odd orders can be shown to 
be zero since the result must be real and r is always ac- 



companied by a factor of i. The result is more compact 
when put in terms of c nm as defined by Eq. I|43|) . In this 
form the result is 

v e (t) = 1-D 4 t 4 -D 6 t 6 + 0(t 8 ), (44) 
Di = 4c 2 nm b 4 nm , (45) 

--De. = c 4 b 6 +c 2 b 6 (46) 

+ CnmC n kbn m bnk \pnk (bnm + &rafc) + ^mfe] , 

where distinct summation over indices is implied for each 
term. Written in this form, we can surmise that the 
r expansion gives convergence for max(6 nm )r < 1 in 
the regime in which most pairs satisfy c nm <C 1. As 
indicated above, typical spectral diffusion problems have 
many pairs that satisfy c nm 3> 1 instead. 

This program has computed algebraic expressions up 
to 0(t 10 ) for both I = 1/2 and I = 3/2. In the 1=1/2 
case, there are 111 eighth order terms, and 1200 tenth 
order terms. Terms of order r 8 have up to four distinct 
nuclear index labels; terms of order r 10 have up to five 
distinct nuclear index labels. The trend that the num- 
ber of distinct nuclear indices is half the order of r will 
be used in Sec. IV B"2l to relate the tau expansion to the 
cluster expansion. 

After computing these algebraic expressions, there is 
an additional challenge of efficiently computing the co- 
efficients of t from these expressions for specific values 
of A n and b nm in order to obtain results such as the r- 
expansion plots of Fig. Appendix \B\ describes how 
this is done. 



2. Relating the tau-expansion to clusters 

The trend seen above is that the maximum number of 
distinct nuclear index labels is half of the order of r. This 
provides a connection to the cluster expansion. By defini- 
tion, a cluster contribution must involve all of the nuclei 
in that cluster. The number of distinct nuclear indices 
for a cluster contribution is thus equal to the cluster size. 
Therefore, clusters of size 2 are, at a minimum, fourth 
order in r; clusters of size 3 are C(r 6 ), etc. Even if we 
cannot prove that this trend continues beyond C(r 10 ) 
we know that clusters of size greater or equal to 5 are at 
least 0(r 10 ). This is a key formal connection between 
the r-expansion and the cluster expansion. 

We can therefore relate orders of A used in Sec. HVI to 
orders of r (up to 5-clusters at least). As long as the r 
expansion is convergent for small N, the cluster expan- 
sion can extend its convergence for large N. As discussed 
above, this convergence in t occurs for max (6 nm )r <C 1 
in the regime in which most pairs satisfy c nm -C 1. We 
also observed, from the argument that ve(t) must be 
real, that no odd orders of r are present. Therefore, as 
discussed in Sec. IIVI the next order of the cluster ex- 
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pansion beyond pairs should include 3-clusters and 4- 
clusters. 



C. Dipolar perturbation expansion 

An expansion in r is not well-suited, on its own, to 
practical problems of interest due to the many nuclear 
pairs satisfying c nm » 1 [defined in Eq. . However, 
in the regime in which c nm 3> 1, a nondegenerate pertur- 
bation expansion in orders of b nm is applicable. We intro- 
duce a bookkeeping parameter A (which we will later re- 
late to the A used in Sec. lIVI making a formal connection 
to the cluster expansion) such that ±Ti ± = Ti° ± XH 1 . 
Here the unperturbed Hamiltonian, Ti° = \ ~}2, n A n I nz: 
is diagonal in the nuclear spin z-basis, while Ti' ~ jTi B is 
the dipolar interaction rescaled to have the same magni- 
tude as Ti°. We note that A here is a purely formal artifice 
to define a perturbation expansion using standard pertur- 
bation theory in quantum mechanics (and later adapted 
to the cluster expansion of Sec. IIVI in order to ensure 
convergence for large systems). 



1. Perturbed Hamiltonian 

Using standard perturbation theory in quantum me- 
chanics, we have the following recursive definitions for 
the eigenvectors and eigenvalues of Ti : 



Plugging this into Eq. (|19f) we have 



ifc±) = \k Q )±\j2\i°) x 

(l°\Ti'\k±) - (l° \k±) (k°\Ti'\k ± ) 



(47) 



E (o) _ E (o) 



K = 



4 0) ±\(k°\H'\k±), 



with the normalization convention that 
(fc°| k) = (k°\ fc°) = 1. 

Now we have 

\k±) (fed 



v ( fc ± l fc ±> 



(48) 
(49) 



(50) 



(51) 



ve(t) 
Cijki 

Fijkl 
Dijkl 

W ijkl 



J2 %L e x p - fc -«« T 



i,j,k,l 
Fijkl 



D 



ijkl 



(I- (i+ \J-) (.1- \K) (K \l-) 
(H K+) (j- li-) (*+ l^+> (I- 



J ijkl 



ijkl ) 



El 



(i°\Ti'\i + )-(j°\Ti' 

- (k°\Ti' \k+) + (l°\Ti' \l-) 



(53) 

(54) 

(55) 
(56) 

(57) 

(58) 

(59) 



We note that C^m and oj'^ kl can be expanded in orders 



of A 



2. Relating the dipolar perturbation expansion to clusters 

In order to show that the dipolar perturbation expan- 
sion can be used as an underlying perturbation theory for 
the cluster expansion (i.e., to supply the A of Sec. HVf) . we 
need to show that for a given term in the A expansion of 
Eq. (|53|l . there are at least as many orders of A as there 
are nuclei involved (cluster size) . In Appendix^] we show 
this to be effectively the case in the max(6„ m )r <C 1 
limit. 

One further point is that ve{t) is an even function of 
A. This is seen by noting that under the transformation 



A 



-A we have U± «-> lit and therefore 



ve{t) = h 
i 



(60) 



{u-u+ulul} 

{uIuIu+U-}\=v b (t). (61) 



This is true because there is a symmetry between "up" 
and "down"; that is, U- <-» L/+ is symmetric within the 
trace operation because I nz — > —I nz V n is symmetric 
within the trace operation. A consideration of finite nu- 
clear temperatures would break this symmetry, but in our 
infinite nuclear temperature approximation there will be 
no odd orders of A. Therefore, as discussed in Sec. lIVI the 
next order of the cluster expansion beyond pairs should 
include 3-clusters and 4-clusters. This was shown to be 
necessary in Sec. lVB"2l when relating cluster sizes to or- 
ders of r, and it is also true here when relating cluster 
sizes to orders of the dipolar perturbation expansion. 



where the summation is over all possible states. From yi CALCULATIONS AND RESULTS 

Eq. U3), 

The information (i.e., the input) necessary to calculate 
U±(t) = - c zfiE ^ T . (52) tne spectral diffusion Hahn echo decay, given in Eq. 1(15)1 

(k± | k±) an( j approximated in the cluster expansion by Eq. I|41|l . 
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are the A n , b nm , and /„ (nuclear spin magnitude) val- 
ues. The internuclear dipolar coupling for b nrn is given 
by Eq. ©. For A n we use 9 - 



An = YlSflWRnW 

1 - 3 cos 2 



(62) 



-e(|i?„ 



ro) 



The first term of Eq. I|62l) is the hyperfine coupling of 
the n th nucleus to the electron, while the second term 
is its dipolar coupling to the electron with r being the 
effective range of the wave function. R n is the position of 
the n th nucleus relative to the center of the electron wave 
function and 9 n is its angle relative to the magnetic field 
direction. While both of these terms were included in our 
calculations, the dipolar contribution of A n proves to be 
negligible for the semiconductor spin quantum computer 
architectures. 



A. Phosphorus donor in silicon 

Our first application is to consider the Hahn echo de- 
cay of the electron spin of a phosphorus donor in natu- 
ral siliconiSi2iii Here is the Kohn-Luttinger wave 
function of a phosphorus donor impurity in silicon, as 
described in Ref. @. Using this in Eq. H62|) , we have 9 - 



A n — 



^ 7s7 f frry \F x {R n ) cos {k Q X n 



(63) 



KFaORiO cos (k Y n ) + F 5 (R n ) cos (k Z n )Y 



" W ^0(11^1 -na), 



Rn 



*i,2(r) 



exp 



+ 



y z +z z 



(na) 2 



\J-K(na) 2 (nb) 



(64) 



with 7s = 1.76 x 10 7 (s G)-\ 7 f = 5.31 x 10 3 (s G) _1 , 
n = 0.81, a = 25.09 A, b = 14.43 A, rj = 186, k = 
(0.85)27r/asi, and asi — 5.43 A. The Si nuclei are located 
on a diamond latticed The central 31 P nuclear spin does 
not contribute to spectral diffusion because its hyperfine 
energy is significantly larger than any of its neighbors, 
suppressing the spin flips by energy conservation. 

In a natural sample of silicon, only a small fraction / = 
4.67% of lattice sites have nonzero nuclear spin. These 
are the spin-1/2 29 Si isotopes, therefore I n — 1/2 for 
all contributing nuclei. We will use (Efc(r)) and (E£(r)) 
to denote Efc(r) and E£(r) averaged, respectively, over 
isotopic configurations with a fraction, /, of 29 Si. We will 
also use the convention that E/.(t) and E|(t) without 
these angle brackets gives the / = 100% result. Thus 



(E*(t)) = / fe S fc (r), 
(E£(t)) = / fe E*(r), 



(65) 
(66) 



where E fe (r), E|(r), and Ej^r) are given by Eqs. (|57Jl . 
(|38H . and 139|) . respectively, taking all nuclei to be 29 Si. 
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FIG. 8: Hahn echo decay ve(t,9) of a phosphorus donor 
electron spin in silicon due to the dipolar nuclear spin bath 
dynamics, (a) Theory (solid lines) and experiment (Ref. IToT) 
is shown for several orientation angles of the magnetic field 
with respect to the crystal lattice, ranging from the [100] to 
the [110] direction (9 = 0, 10, 20, . . . , 90). (b) Here we plot 
— In ve (t, 0) + In ve (t, 8 — 0) , allowing for the removal of any 
decoherence mechanism which is independent of 6. The quali- 
tative and quantitative agreement between theory and exper- 
iment is remarkable, in contrast to the stochastic approach 
(dashed). 



The fact that only a fraction, /, of these nuclei con- 
tribute to the diffusion is accounted for by the f k factors 
in Eqs. I|65(l and (|66|l because f k is the probability that 
all nuclei in a cluster of size k have nonzero spin. 

For the spin-1/2 nuclei that contribute to spectral dif- 
fusion in natural silicon, we can write the following ana- 
lytical solution for pairs (2-clusters) using Eq. 119|) : 



(i 



[cos (oj nm r) - lY 



2b nm y/l + C 



2 

run ■ 



(67) 
(68) 



where c nm is given by Eq. (|43ll . 

Our numerical calculations of Hahn echo decay in 



the lowest order of the cluster expansion, 



(2) 



(r) 



exp((E 2 (r))) using Eqs. (JSSJl, (123, and JBJJl, are 

shown for several magnetic field orientation angles in 
Fig. Eta) with a direct quantitative comparison to the 
experiment. 19 The dipolar coupling [Eq. @] contains 
an important anisotropy with respect to the 9 nm an- 
gle formed between the applied magnetic field and the 
bond vector linking the two spins (_R nm ). This property 
leads to a strong dependence of spin echo decay when 
the sample is rotated with respect to the applied B field 
direction. The experimental data is taken for bulk natu- 
ral silicon with phosphorus doping concentration equal to 
2 x 10 15 cm 3 |10|. The high concentration of phosphorus 
donors leads to an additional decoherence channel arising 
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from the direct spin-spin coupling between the electron 
spins that contribute to the echo. This contribution can 
be shown to add a multiplicative factor exp (— r/1 ms) to 
the Hahn echo^ Because this contribution is indepen- 
dent of the orientation angle, we can factor it out by sub- 
tracting the 9 — contribution from the logarithm of the 
experimental data taken at angle 9. The result is shown 
in Fig. Etb) (log-log scale). Our theory seems to explain 
the time dependence of the experimentally observed echo 
quite well. This result is to be compared with the recent 
stochastic theory of Ref.0 [Dashed line in Fig.[8Jb) shows 
the stochastic calculation for 9 = 60°]. Although the 
stochastic theory yields roughly correct coherence times 
in order of magnitude, it fails qualitatively in explain- 
ing the time dependence [that is, the shape of the decay 
as can be seen from the incorrect slope of the stochastic 
calculation in the log-log plot of Fig. lEfb)]. The present 
method is able to incorporate all these features within a 
fully microscopic framework, obtaining both qualitative 
and quantitative agreement with experiment. 

An important issue in the context of quantum infor- 
mation processing is the behavior of spin coherence at 
the shortest time scales. The experimental dataiS in 
Fig. [S] reveals several oscillatory features which are not 
explained by our current method. These are echo mod- 
ulations arising from the anisotropic hyperfine coupling 
omitted in Eq. (JSJ). 11 This effect can be substantially re- 
duced by going to higher magnetic fields. (In a quantum 
computer B ~ 9 Tesla will probably be required in order 
to avoid loss of fidelity due to echo modulation^ 4 *) On the 
other hand, spectral diffusion itself is essentially indepen- 
dent of magnetic field even to extremely high magnetic 
field values (B ~ 10 Tesla). The echo modulation effect 
is expected to be absent in III-V materials^ hence our 
theory allows the study of spin coherence at time scales 
of great importance for quantum information purposes, 
i.e., very short time scales, as long as echo modulation 
effects are quantitatively unimportant. 

Isotopic purification can reduce the value of / (fraction 
of 29 Si nuclei). Figure [5] contains information that is use- 
ful for understanding how the Hahn echo curves change 
as / is changed (i.e., lowered via isotopic purification). 
In a log-log plot, In (ve{t)) ~ (E 2 (t)) oc / 2 simply shifts 
vertically when / is changed. Figure [5] shows both the 
/-independent E 2 (r) (i.e., / = 100%), and (E 2 (r)) for 
natural Si (/ ~ 5%). Results are shown for magnetic 
field angles that yield the extremal slowest and fastest 
decoherence. For natural Si, in a wide range of r about 
Ti/ e , where Ve(tiu) — 1/e, (E 2 (t)) matches r 23 curves 
very well. In this range of r, therefore, we may write 

v E {r) « exp(-/ 2 (r/r ) 2 - 3 ) (69) 

= cxp(-(r/r 1/e ) 2 - 3 ), (70) 

where 

T\/e = TO// ' OC / , (71) 

providing a formula that allows us to adjust our Hahn 
echo curves to other values of / for a range of r in which 
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FIG. 9: Lowest order results for the natural log of the Hahn 
echo, ln(ti_B(r)) « (E2(r)) oc f 2 , for Si:P in a log-log plot. 
The solid lines give £2(1") with / = 1. Dotted lines give 
<E 2 (t)> for natural Si (/ = 4.67%). In this log-log plot, mul- 
tiplying by f 2 simply shifts the curves vertically. Isotopic pu- 
rifucation would shift these curves up further. The two mag- 
netic field angles shown give extremal results. Corresponding 
to 6 angles in Fig. B || [100] is 9 = 0° and B || [111] is 
8 « 54.7°. Dashed lines fit the natural Si curves near their 
— 1 values (where ve ~ 1/e) with t 2 ' 3 power law curves (linear 
in the log-log plot). 

Eq. Ij69(l is applicable. Shortly after the original submis- 
sion of this manuscript, experimental data from Abe et 
al. appeared in the literature 3 ^ which shows echo time 
scaling ranging between / _0 - 86 and / _0 ' 89 , in remarkable 
agreement with our prediction [Eq. I[71jl]. Tyryshkin et 
alw& report Si:P Hahn echo decay forms of exp (r 2 ' 4 ), 
in agreement with Eq. I|69(l . with exception to magnetic 
field orientations near the [100] direction. In the [100] 
direction, they report a form of exp ( r 3 ± - 2 ); the reason 
for this discrepancy remains unknown. 

We now check the convergence of our cluster expansion 
for this Si:P system. Using Eq. Ij41(l and averaging over 
the isotopic configurations yields 

k 

This approximates the ideal cluster expansion [see Sees. 
IIV B 21 and IIV B 3j with an error that we may estimate 
as (E*(r)) = (E|(r)) + (E|(r)). This error is esti- 
mated by the correction needed to compensate for over- 
lapping pairs [either the same pair overlapping itself, 
E 2 (t), or two different pairs overlapping, Eg(r)] in the 
approximation. Figure 1101 shows these relative correc- 
tions, E^(r)/E 2 (r) and E|(t)/E 2 (t), to \u(v e (t)) for 
both / = 100% and natural Si (/ - 5%). The graphs 
also show the respective Hahn echoes (in a log time scale 
which alters their appearance) to show that these relative 
corrections are very small up to the tail of their respective 
echo decays. The argument, given in Sec. IIV B 31 that 
E* (t) would be small was only applicable for Ve(t) > e _1 
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FIG. 10: Relative errors (with scales on the right) to the 
log of the Hahn echo due to overlapping pairs for both 100% 
29 Si (top graph) and natural Si (bottom graph). In these ex- 
amples, B || [100]. The Hahn echoes themselves are shown, 
as well, with the to 1 scales on the left. All curves share 
the same logarithmic time (r) scale. It is apparent that these 
relative corrections are very small up to the tail of their re- 
spective echo decays. 



so it is expected that this approximation approaches fail- 
ure out in the tail of the decay. This is irrelevant for 
practical purposes. 

The expansion of Eq. (|72|l is convergent where 
(£ fe+1 (r)) < (Sfe(r)) (implying that XL < 1 effec- 
tively). The Sfc(r) functions have been calculated (up 
to k = 5) using Monte Carlo techniques with cluster con- 
tributions, v' c (t), for clusters that are larger than pairs, 
calculated by numerically diagonalizing H ± [Eq. I|18l) ]. 
For each (t) independently, the maximum distance be- 
tween neighbors and the maximum distance of nuclei to 
the donor is increased for various Monte Carlo runs un- 
til convergence within a desired precision is reached. To 
speed up each Monte Carlo run, clusters are chosen with 
a heuristic bias for those that have strong coupling be- 
tween the constituent nuclei as well as a bias for clusters 
closer to the donor. Appropriate weighting factors are 
used to counteract these biases. 

Figure ITT1 compares /-independent (i.e., / = 100%) 
Efc (t) functions in a dual (showing positive and negative 
values) log-log plot for Si:P with B || [100]. In other 
words, it compares successive orders of the expansion for 
the natural log of the Hahn echo, \ii(ve(t)), with the 
/ dependence removed. As one might anticipate by the 
fact that £ 3 (t) and S 4 (r) are both 0(A 4 ) (Sec.|V|proved 
that there are no odd orders of A for either perturbation 
theory), they are similar orders of magnitude, at least for 
the 0.03 ms < t < 1 ms range. Near r ~ 1 ms, however, 
the perturbation theory fails [having the condition that 
max(& IlTO )r <C 1] as we see that |E4(r)| surpasses |E 3 (r)|. 
Interestingly, all orders approach the same order of mag- 
nitude near r ~ 1 ms. We can thus identify the break- 
down of the cluster expansion. Note, however, that this 
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FIG. 11: Successive orders of the cluster expansion for the 
natural log of the Hahn echo [Eq. 1721 1. computed for Si:P 
and B | [100], for the 100% 29 Si theoretical scenario. The 
thick black line gives the lowest order result, £2(1"), and other 
solid lines give higher order Sfc(r) results. The dotted lines 
give the negative of their corresponding functions provided to 
assist in the absolute value comparison of these higher order 
corrections. A failure of convergence occurs near r ~ 1 ms 
where all of the curves are the same order of magnitude. This 
occurs well into the tail of the decay, however, and therefore 
has no practical consequence. 



is well into the tail of the decay [where ve{t) < e~ 1000 ] 
and therefore this breakdown is irrelevant for practical 
purposes. It is prudent, in any case, to understand the 
limitations of this expansion. 

The (S^(t)) curves for some fraction, /, of 29 Si will 
be the same as the £fc(r) curves in the log-log plots of 
Fig. II II except with appropriate vertical shifts [multiply- 
ing by f k effectively appears as addition by k log (/) in 
the log plot] due to the / dependence. Higher orders will 
be shifted closer to zero than the lower order curves and 
therefore these curves will be more separated (actually 
improving the cluster expansion convergence). The top 
graph of Fig. ^| is analogous to the bottom (negative 
range) graph of Fig. ITTl for natural Si (dashed lines indi- 
cate negated curves, i.e., where values are actually posi- 
tive) . We show only the low order corrections to the log of 
the Hahn echo, including (E*(r)) as well as (E 3 (t)} and 
(£4(7-)) and not bothering with (E 5 (t)). (E 3 (t)), with 
its inclusion of 3-cluster, gives the largest correction. Al- 
though E*(t) is of a comparable order of magnitude, its 
correction partially cancels the £ 3 (r) correction because 
they are opposite in sign. We may therefore use E 3 (r) for 
a conservative estimate of the error of the lowest order 
cluster expansion result. The bottom graph of Fig. IT51 
shows the absolute (as opposed to relative) error of the 
lowest order Hahn echo result estimated by the inclusion 
of E 3 (t). The Hahn echo is displayed for reference. At its 
maximum, this absolute error is approximately 0.001. Al- 
though our cluster expansion fails near r ~ 5 ms [where 
|(S 2 (r))| ~ |<E*(r))| ~ |(E 3 (r))| ~ |(E 4 (r))|], the abso- 
lute error will stay small if we assume that our Hahn echo 
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FIG. 12: Low order corrections to the Hahn echo, ve(t), com- 
puted for Si:P in natural Si (/ = 0.0467%) with B || [100]. 
(top) Log-log plot of low order contributions to the natural 
log of the Hahn echo, In (ve(t)). Ordinate axis is negative as 
in the bottom graph of Fig. HTI however, dashed lines indicate 
negated curves (and thus represent positive values), (bottom) 
Conservative estimate of the absolute error of the lowest order 
Hahn echo result (scale on the right) due to 3-cluster contri- 
butions, (£3(7-)). The lowest order Hahn echo result is shown 
as a reference (scale on the left). The logarithmic time scale 
is the same for all plots (top and bottom). 
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FIG. 13: This compares the r-expansion results with the low- 
est order cluster-expansion result, exp[£2(r)] [Eq. I|320 ]. for 
an arbitrary set of nuclei (ranging from 15 to 20 lattice con- 
stants away from the center) for the phosphorus donor in nat- 
ural silicon. Shown are the r-expansion results at the highest 
two orders of r that we have computed. Note the agreement 
with the cluster expansion results before the two T-expansion 
results diverge. Also shown is 1 + £2(7"); its disagreement in- 
dicates the importance of accounting for simultaneous pairs 
in the cluster expansion. 



decay is forever monotonically decreasing. For all practi- 
cal purposes, the lowest order result is therefore valid up 
to 0.1% of the initial Ve(0) = 1, and higher order terms 
only provide corrections beyond 99.9% accuracy level. 

We have also verified that our cluster expansion results 
agree quantitatively with the r-expansion [Eq. 144|) ] (ap- 



plied directly, not in the context of the cluster expansion) 
up to tenth order in r when excluding nuclei close to the 
center of the electron wave function where c nm 1. 
This provides an independent verification of the cluster 
method approach. An example is given in Fig. 1131 As 
discussed in Sec. IV Al however, numerical calculations 
show that dipolar perturbation theory accounts for the 
Hahn echo decay in physically relevant systems we have 
studied. This is because clusters for which this perturba- 
tion theory is applicable will tend to operate on shorter 
time scales (due to the corresponding large hyperfine en- 
ergies) . This has been shown by comparing the lowest or- 
der cluster expansion using exact pair contributions ver- 
sus approximate pair contributions in the lowest dipolar 
perturbation order as discussed in Sec. IV Al 



B. Gallium arsenide quantum dot 

Next we consider the Hahn echo decay of a localized 
quantum dot electron spin in GaAs. For this, ^(Rn) 
will be the quantum dot wave function parameterized by 
the quantum well thickness, z , and Fock-Darwin radius, 
£{B) (a function of the magnetic field strength), as de- 
scribed in Ref . @. Using this in Eq. (|62|l , we have^ 



16 7 s7/%j As /4) 
3 £ 2 {B)z a 

Xl + V? 



d(I)cos^ —Z n 



(73) 



x exp 



£ 2 (B) 
1 - 3 cos 2 6 n 



&{z /2-\Z n \) 
@[Xl + YZ-e{B)], 



with a GaAB = 5.65 A and 75 = 1.76 x 10 7 (s G)" 1 (the 
free electron gyromagnetic ratio). The GaAs lattice has 
a zinc-blende structure with two isotopes of Ga atoms 
placed on one fee lattice and 75 As atoms placed on the 
other fee latticed The Ga isotopes are 60.4% 69 Ga 
and 30.2% n G&2£L We used 7/ = 4.58,8.16,6.42 x 
10 3 (s G)- 1 and d{I) = 9.8, 5.8, 5.8 x 10 25 cm" 3 for 75 As, 
71 Ga, and 69 Ga, respectively All of these nuclei have 
a valence spin of / = 3/2 which means that Eq. (|67|l is not 
applicable. Instead, cluster contributions are computed 
by numerically diagonalizing 7i . 

Most of our results only include dipolar interac- 
tions for b nm [Eq. ©]. However, indirect exchange 
interactionaSSiSL 2 ^ between the nuclei of GaAs may be of 
the same order of magnitude as the dipolar interactions 
[Eq. @] between nearest neighbors. 13 There is enough 
quantitative ambiguity in the literature about the indi- 
rect exchange interaction that we have chosen to leave 
it out of our calculations in order to have a precise the- 
ory for the dipolar nuclear spin bath dynamics only. We 
make an exception near the end of this section where we 
make a comparison to the results in Ref. ^3 that include 
this interaction. 

With exception to the I nz I m z term in H B [Eq. ©], 
the different types of nuclei are decoupled. In the low- 
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est order of the cluster expansion, the different types of 
nuclei truly are decoupled because a pair cannot con- 
tribute to the spectral diffusion without the flip-flop in- 
teraction between them. Furthermore, this b nm I nz I mz 
term is negligible where the dipolar perturbation is ap- 
plicable because the diagonal of the Hamiltonian in the 
nuclear z-basis will be dominated by TL A . We will later 
discuss comparisons made between calculations using ex- 
act pair contributions versus approximate pair contribu- 
tions using the lowest order dipolar perturbation (also 
done in Sec. I VI Al for Si:P) to demonstrate that the dipo- 
lar perturbation theory is applicable. Using these argu- 
ments, we will simplify our calculations by approximating 
Eq. © with 



i 'Jm ) Ai+ Im - 



(74) 



This approximation effectively decouples the As nuclei 
from the Ga nuclei and the Ga isotopes from each other. 
Although it is possible to perform our cluster calcula- 
tions without this approximation (we have done so for 
some lowest order calculations with no additional com- 
plication and the difference is indeed negligible), it avoids 
unnecessary complications for higher order calculations. 
We may then perform calculations for each of these nu- 
clear isotopes separately. Using the subscript x to rep- 
resent each of these types of nuclei, the contributions to 
the natural log of the Hahn echo become 
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(75) 
(76) 



where E fex (r), S| x (r), and SJ x (r) [Eqs. f^, (HJ|, and 
(|39|l . respectively] sum over clusters in the appropriate 
sublattice of isotope x, and f x is the fraction of this iso- 
tope in its sublattice. For 75 As, 71 Ga, and 69 Ga, it is 
appropriate to use f x = 100%, 30.2%, and 60.4%, respec- 
tively. Monte Carlo techniques (described in Sec. IVI A(l 
are used to approximate these sums. 

The lowest order results, ve(t) ~ exp [(£ 2 (t))], for 
most of our GaAs calculations show a Hahn echo decay 
of the form exp [— (2r/to) 4 ] • This differs qualitatively 
from the decay for the Si:P which, by our calculations, is 
in the form exp [— (2r/to) Q ] where a ~ 2.3 for a range of 
t appropriate for natural Si and some range of isotopic 
purification. The form of the GaAs echo decay does not 
change if we repeat the calculation with 7 = 1/2 rather 
than 3/2, although to does change. The qualitative dif- 
ference between Si:P and GaAs dot solutions is therefore 
not due to the difference in spin magnitude of the nu- 
clei, but rather the difference in the form of the electron 
wave function and lattice occupation leading to different 
distributions of c nm values. GaAs, in ranges of tested 
zq and £ parameters, has a more narrow distribution of 
Cnm for significantly contributing pairs than that of the 
Si:P problem. As a result of this narrow distribution, the 
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FIG. 14: For GaAs quantum dots, to (circles), 
istic initial decay time, and t\i e (diamonds), the e _1 decay 
time, vs the Fock-Darwin radius I for various quantum well 
thicknesses, zq = 5, 10, and 20 nm. The orientation of the 
magnetic field is (a) along the Zo confinement of the quantum 
dot and [100] lattice direction, or (b) perpindicular to the zo 
confinement direction and along [110] of the lattice. 



fact that the lowest nontrivial order of the r-expansion 
is 0(t 4 ), which is universal for all values of /, expresses 
itself in the sum of pair contributions in £2(1") [Eq. (|37() ]. 
This exp [— (2t/£o) 4 ] form is confirmed by Yao et al*& 
using a completely different theoretical technique. Yao 
et al. use a pair-correlation quasiparticle approach to 
the "spectral diffusion" problem that is equivalent to our 
lowest order cluster expansion, and they give results for 
quantum dots in GaAs. They come to the same conclu- 
sion that the form of the decay depends on the relative 
breadth of the excitation spectrum. 

Figure PHI shows the to of the initial exp [— (2r/io) 4 ] 
Hahn echo decay for various parameter settings of zq 
and £ with two different magnetic field orientations. Also 
shown is ti/ e , defined such that ve{t — ti/ e /2) = e _1 . 
One can think of this t 1 / e as an effective T2-time for the 
problem although the echo decay is not a simple expo- 
nential. Except for small dots, to = ii/ei indicating that 
the decay has the form exp [— (2t/to) 4 ] . Small dots de- 
viate from this form, beginning to have longer t\/ e decay 
times than their initial characteristic times, to- It w &s 
noted in Ref. |9| that decoherence times become infinite 
as the size of the quantum dot approaches zero or infin- 
ity with a minimum decoherence time at some finite size. 
The former is simply because the electron has no inter- 
action with nuclei as the quantum dot size approaches 
zero, and the latter is because the nuclei all have the 
same coupling to the electron as the size becomes infi- 
nite. For zo = 5 nm we begin to approach this maximum 
decoherence (minimum ti/ e ) near £ = 10 nm, but only in 
the regime where t\i e deviates from to- 

As discussed previously, the Ga and the As nuclei, on 
separate fee lattices, are decoupled by our approxima- 
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FIG. 15: For GaAs quantum dots, to, the characteristic initial 
decay time, and ii/ e , the e _1 decay time, vs the Fock-Darwin 
radius £ for a quantum well thicknesses of zo = 2.8 nm. The 
orientation of the magnetic field is parallel to [110] of the lat- 
tice; this is perpendicular to the zo confinement of the quan- 
tum dot. This shows results both when including or excluding 
the indirect exchange coupling. 



tion [Eq. 1)74(1 . In silicon, the asymmetry of the diamond 
lattice results in maximum decoherence in the [111] direc- 
tion. In this case, because the fee lattice is more symmet- 
ric, the angular dependence is primarily a result of the 
shape of the quantum dot (not the lattice). Figure IT"T1 
shows slight quantitative differences when the magnetic 
field is along the Zq confinement direction or perpendic- 
ular to it. 

Because most of our GaAs results are in the form cor- 
responding to the limit of small r, it is tempting to think 
that GaAs is dominated by the c nm <C 1 regime appropri- 
ate for the r-expansion. However, as with the phospho- 
rus donor in Si, we have compared calculations that use 
exact pair contributions versus approximate pair contri- 
butions using the lowest order dipolar perturbation (also 
discussed in Sec. IV AJ) . These different calculations agree 
very well for small quantum dots, but deviate slightly 
for larger quantum dots. Intermediate sized dots are 
well-approximated by either perturbation theory. This 
is merely a result of small decoherence times such that 
the r approximation is valid even though the dominating 
pairs have c nm > 1. The agreement with the perturba- 
tion expansion gives justification for the approximation 
used in Eq. I|74[l in which b nm I nz I mz was neglected. 

We now return to a discussion of the indirect exchange 
interaction between nuclear spins (mediated by virtual 
interband electronic transitions) that were neglected in 
the above calculations. Including the exchange interac- 
tion, we should use 
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FIG. 16: Low order corrections to the Hahn echo, ve(t), com- 
puted for a GaAs quantum dot with B || [110], zo = 10 nm, 
and £ = 50 nm. (top) Log-log plot of low order contribu- 
tions to the natural log of the Hahn echo, ln(ii_B(r)). Ordi- 
nate axis is negative; however, dashed lines indicate negated 
curves (and thus represent positive values), (bottom) Conser- 
vative estimate of the absolute error of the lowest order Hahn 
echo result (scale on the right) due to 3-cluster contributions, 
(£3(7-)). The lowest order Hahn echo result is shown as a 
reference (scale on the left). The logarithmic time scale is the 
same for all plots (top and bottom). 



Si:P system to a high degree of accuracy. Yao et a/J 3 . per- 
formed spectral diffusion decoherence calculations (using 
an equation that is equivalent to our lowest order result) 
for GaAs quantum dots including the indirect exchange 
interaction. As a verification of the correctness of our 
calculations, Fig. reproduces their Hahn echo results 
using our method but including the indirect exchange 
interactions as given bji22 



b e = 

nm 



"7n7« 
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(78) 



where b^ m is the dipolar coupling [Eq. ©], and b^ m is the 
indirect exchange coupling. We note that b e nm = in the 



There is ambiguity in the literature about the appropri- 
ate sign of Eq. (|78fl . We have adopted the convention 
used by Yao et ali& in order to reproduce their results. 
Figure H31 also shows the results for the same parameters 
when the indirect exchange is excluded; it is apparent 
that this coupling is significant, at least in GaAs quan- 
tum dots for these choices of the exchange coupling. The 
kink in the t±/ e curve for the case of excluded indirect ex- 
change is believed to be a discrete lattice effect only no- 
ticeable for small quantum dots. For such small quantum 
dots, it is likely that Eq. J73J), derived from an approxi- 
mate electron wave function, is somewhat inaccurate. 

With our cluster expansion approach, we can estimate 
the error of our calculated decay curves by performing 
higher order calculations. We observe that the larger 
quantum dots have larger corrections. For quantum dots 
with zq = 20 nm and £ — 100 nm our calculations indi- 
cate maximum correction to the Hahn echo decay curves 
on the order of 10 -3 , 0.1% of the initial u_e(0) = 1, just 
as it was for natural silicon (Fig. I12fl . For dots with 
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zq = 5 nm and £ = 10 nm, absolute corrections are on 
the order of 10" 4 , 0.01% of the initial v E (0) = 1. Fig. [HO 
analogous to Fig. El gives these low order corrections ex- 
plicitly for an intermediate size (zq = 10 nm, £ = 50 nm) 
GaAs quantum dot. 



VII. CONCLUSION 

In conclusion, we describe a quantum approach for the 
problem of localized electron spin decoherence due to 
the fluctuation of dipolar coupled nuclear spins. In con- 
trast to former theories, our method requires no ad hoc 
stochastic assumption on the complex dynamics of the 
environment responsible for decoherence. Hence it pro- 
vides an important example where direct integration of 
the environmental equations of motion provides a system- 
atic understanding of the loss of coherence which needs 
to be controlled for quantum information applications. 

The most important theoretical accomplishment of our 
work is the development of the first fully quantum mi- 
croscopic theory for the localized electron spin decoher- 
ence due to the spectral diffusion induced by nuclear spin 
bath dynamics. The important nuclear spin dynamics 
for the spectral diffusion problem is the dipolar inter- 
action induced nuclear spin flip-flops although for the 
GaAs quantum dot system, we have also included the 
indirect exchange interaction between the nuclear spins 
which turns out to be quantitatively comparable to the 
dipolar contribution to spectral diffusion. Our results are 
formally exact, and our numerical calculations provide 
an essentially exact (to better than 0.1% of the initial 
value) quantitative description of the Hahn spin echo de- 
cay. The significance of our quantum theory lies in the 
fact that, unlike all other theoretical descriptions of spec- 
tral diffusion spanning the last 50 years, we do not make 
any ad hoc phenomenological stochastic approximation 
in dealing with the non-Markovian spin dynamics in the 
spectral diffusion phenomena. We solve the problem es- 
sentially exactly using a quantum cluster decomposition 
technique, which is then theoretically justified by car- 
rying out calculations to higher orders and by compar- 
ing with two independent perturbation techniques (i.e., 
the tau expansion and the dipolar perturbation theory). 
Very recently, a completely independent verification of 
our theory and results has appeared in the literature M> 

Our theory for spin decoherence is restricted to under- 
standing the so-called T 2 -decoherence of a single localized 
spin in a solid due to the dipolar nuclear spin flip-flops 
in the surrounding nuclear spin bath and in the specific 
context of a spin echo experiment. We neglect spin de- 
coherence effects arising from the direct hyperfine cou- 
pling between the central spin and the surrounding nu- 
clear spins, which would contribute to the free induction 
decay signal, but can be corrected by the spin echo tech- 
nique. Our method uses a cluster expansion technique 
supported by two underlying perturbation theories: the 
dipolar perturbation expansion and the tau-expansion. 



Each of these perturbation theories generates a unique 
expansion scheme in its own right, but only the clus- 
ter expansion converges for large systems [only for short 
times, r -C max(6 nm ) , but often extending to the tail 
of the decoherence decay]. The convergence of this ex- 
pansion has to be explicitly checked in each application, 
as we have done in this work. It is possible that an al- 
ternative more conventional perturbation theory for the 
non-Markovian spectral diffusion problem can be devel- 
oped within the standard interaction picture of the many- 
body quantum theory which would be completely equiva- 
lent (for short times) to our cluster expansion technique. 
More work would be needed to clarify and establish such 
a possibility. 

We note that spectral diffusion induced decoherence 
cannot be characterized by a simple T2 time since the 
Hahn echo decay is not a simple exponential, and in 
fact, obeys different temporal power laws in the Si:P and 
the GaAs quantum dot systems. Spectral diffusion is a 
pure dephasing process arising from the temporal (non- 
Markovian) magnetic field fluctuations at the localized 
electron spin due to the nuclear spin dynamics, but it 
cannot be parameterized by a single T2 time except as a 
crude approximation. Within this crude approximation, 
we find the spectral diffusion induced Ti to be about 
100 /its for the Si:P system and about 10 /is for the GaAs 
quantum dot system. But, this T2 can be enhanced in- 
definitely (up to tens of milliseconds) in the Si:P system 
through the isotopic purification of Si (i.e., by removing 
29 Si nuclei from the system) whereas in the GaAs quan- 
tum dots, T2 ~ 10 /xs is essentially an absolute upper 
limit (when using simple Hahn echo refocusing) since all 
Ga and As nuclei isotopes have free spins contributing to 
the spectral diffusion and isotopic purification is impossi- 
ble. It is important to emphasize here that although spin 
polarizing the nuclei (e.g., through the dynamic nuclear 
polarization technique) would, in principle, suppress the 
nuclear induced spectral diffusion decoherence of electron 
spin, in practice, this would lead only to rather small en- 
hancement of electron spin coherence since the presence 
of even a few nuclei with the "wrong" spin would cause 
nuclear pair flip-flop processes^ 

Finally, we comment on the fact that the spectral dif- 
fusion process is quite a generic and general phenomenon 
in any spin decoherence problem with coupled spin dy- 
namics (e.g. electron and nuclear spins, different types 
of nuclear spins, etc.) where the dynamics of one spin 
species has nontrivial (i.e., non-Markovian) temporal ef- 
fects on the evolution of the spin dynamics of the other 
species. For example, a trivial (but not often empha- 
sized in the literature) consequence of spectral diffusion 
consideration implies that in systems (e.g., Si:P; GaAs 
quantum dots) of interest to quantum computer archi- 
tectures, the single flip of the localized electron spin will 
immediately decohere all the nuclear spins in its vicinity. 
Thus the nuclear spin Ti time in these systems can at 
most be the T\ time for the electron spin! The typical 
low-temperature T\ time for electron spins in the GaAs 
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quantum dots has been measured to be 1 ms or so, and 
therefore the nucleus spin T 2 time would at most be 1 ms 
in the GaAs quantum dots, at least in the neighborhood 
of the localized electrons in the dot. The same consid- 
eration applies to the Si:P system. We believe that the 
general quantum theoretical techniques developed in this 
paper will be helpful in the studies of the temporal dy- 
namics of other coupled spin systems wherever one spin 
species could act as a "decoherence bath" for the other 
system. 
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APPENDIX A: HAHN ECHO ANALYSIS 

Here we derive Eq. I|16|l from Eq. Hll|) and related equa- 
tions. By noting that 



U (t) = e-^e-^U^t), 



-i{H A +H B )t 



Using 



{cr x , e , S z } = =>■ a x , e e 



(A5) 
(A6) 

(A7) 



and the fact that H. Zn commutes with a X:S (the former 
acting only on the nuclear spins and the latter acting on 
the electron spin), we can rewrite Eq. 113|) 



U(t) = U {T)a x>e U {r) 

U'(t) = U' Q {r)a x Mr). 
We can then rewrite Eq. (|T2|l 



in z "(2r) 



(A8) 

(A9) 

(A10) 
(All) 



P (t) = U'(r)e- mZn ^p a e mZn ^U'^T) (A12) 
= U'(t) P0 U>Ht), (A13) 

using the fact that Tt Zn commutes with po (since it com- 
mutes with \x° e ) and H n = H z » + U B ). 

At this point it is convenient, for clarity, to write our 
operators as block matrices in the electron spin z basis. 
First we have, from Eq. (IA6I) 



U+ 
U- 



where 



S x + iS y = S + = \l e )(O e \, (Al) 

we may rewrite Eq. 111|) . 

v b (t) = 2 Tk n «O e |p(r)|l e )), (A2) with 

where the n subscript in Tr„ denotes a trace over the 
nuclear subspace. In addition to this equation, we will be 
referring to Eqs. CEJ-©, ©i an d JHl f° r the free evolution 
Hamiltonian, H = H Ze + H Zn + H A + H B , as well as 
Eqs. I|12|) - (|15fl t° define the density matrix, p(r), in terms 
of the evolution operator, U(t), and inital density matrix, 
Po(r). 

First we note that Tt Ze commutes with the rest of the 
Hamiltonian and that TL Zn commutes with TL A (since an 
operator commutes with itself and any operator that acts 
on an orthogonal state space) . TL Zn also commutes with 
7i B (and therefore commutes with the entire Hamilto- 
nian). This is proven by noting that Ti. B preserves the 
overall polarization of like nuclear spins. That is, 

[jnlnz + Imlmz, Kin ~ 7m)4+4-] = 0. (A3) 
Let us define 

U (t)=e- mt , (A4) 
and exploit the above commutation relations to write 



U ± = e- m r ; U\ = e m ±\ 



n B ±n An , 

-^A I 



(A14) 

(A15) 

(A16) 
(A17) 



U± can be interpreted as the evolution operators of the 
nuclei with the electron spin constrained to be up or 
down, respectively, and ignoring the external magnetic 
field. Although U+ and [/_ are functions of r, we have 
dropped this as an explicit parameter for convenience. 
In this same electron spin basis, we can write 
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(A18) 



in order to obtain, from Eqs. (fl^l . 1A13|1 . (|Alip . and 

EH , 



(0 e |p(r)|l e 



e i4 
2M ' 



:U-U + e- nn/T U f _Ul. 



(A19) 



Plugging this into Eq. I|A2|) gives 



ve(t) 



e i<t> 
~M 



Tr {U-U+e~ H "/ kBT U f _Uiy (A20) 



We are interested in the echo envelope which is given by 
the magnitude of Eq. (|A20|) and thus the e 1 ^ phase factor 
may be dropped. 



22 



APPENDIX B: INDEPENDENCE OF CLUSTER 
CONTRIBUTIONS 

In order to prove the unique existence of the decom- 
position of vs(t) [the solution to the Hahn echo, i>_e(t), 
when only including the nuclei in some set S] given by 
Eq. (|20|) . we must prove that each cluster contribution 
is independent of anything outside of that cluster. Only 
then can we assume that the contribution of a given clus- 
ter is the same in every instance (i.e., in every term that 
involves this cluster) in such a decomposition. 

Let us consider some cluster C contained in S. Because 
we are only concerned with the independence of processes 
in disjoint sets, let us remove any processes from vs(t) 
with coupling between nuclei in C and nuclei in S — C. 
Define v_a-.b( t ) to be the sum of all of the terms in the 
sum of products expansion (described at the beginning 
of Sec. lIVB 1(1 of V(A[j 8) i T ) containing any factor of b nm 
in which n £ A and m £ B (or vice versa). Thus, by 
definition, vs(r) — vc-.(s-C)( T ) does not involve any cou- 
pling between a nucleus in C and a nucleus in S — C; it 
only involves coupling amongst nuclei in C or amongst 
nuclei in S — C. This is the same as calculating Vs(t) if 
we artificially impose the condition that 



b nm = b mn — V n £ C, m £ S -C. 
For convenience, wc define the following: 

V s = u s -u s+ u s _ul + , 



U s ± = e 



-VHtr 



1 



n£S 



such that 



v S (t) = ± Tr {V S } = Tr 5 {V s }, 



(Bl) 

(B2) 
(B3) 
(B4) 



(B5) 



where Ms is the number of composite states for the nuclei 
in set S and Tr,s takes the trace over only these states. 
Note that if we artificially impose condition I|B1|) , then 



TL C + TL 



Noting that Ti^ commutes with 7~tf s _c) (since they act 
on disjoint sets of nuclei), then 



± 

(s-cy 
± 

(S-C) 



(B6) 



U s ± -» e- m c T e~ in is-c) T (B7) 

= U c ±U { s-c)±, (B8) 

V s - VcV ( s-c)- (B9) 

We may therefore write 

vs(t) - v C: (s-c)(t) = tt Tt s {VcV { s-q} (BIO) 

Ms 



Tr c {T/ c }Tr 



(S-C) 



{%-c)| 



Bll) 



M c M (S -c) 
vc(t)v {S -c){t). (B12) 



Equation IjBlip is a direct consequence of the fact that 
Vc and V($-c) ac t on orthogonal subspaces of the Hilbert 
space because they act on disjoint sets of nuclei. We 
conclude that 



Vs{t) = V C (t)v (S -C)(t) +V C :(S-C)(t), 



(B13) 



proving that processes involving disjoint sets of nuclei 
are independent for the following reason. Any interde- 
pendent process of cluster C must be contained in the 
vc (t) part of the right-hand side of Eq. (|B13|) . This has 
no dependence on the set S beyond the requirement that 
C be contained in it. For this reason, the cluster contri- 
bution for cluster C is well-defined in the sense that it is 
independent of all other nuclei. 

We may also verify from Eq. (|B13|) that we are to mul- 
tiply independent cluster contributions in order to ac- 
count for simultaneous processes [the first term on the 
right side of Eq. (|B13p ] and that alternative combina- 
tions are to be included via addition (the second term). 
With this observation, we can decompose any vs(t) into 
a sum of products of cluster contributions. Recall that 
a cluster contribution includes all contributions from in- 
terdependent processes involving only and all of the nu- 
clei in the cluster. Each process is thus assigned to only 
one cluster and there is no redundancy or overlap be- 
tween different cluster contributions. Therefore, without 
redundancy, vs(t) is the sum of all possible products of 
cluster contributions as given by Eq. lf2*U|) . 



APPENDIX C: METHODS FOR COMPUTING 
THE ALGEBRAIC TAU-EXPANSION 

We describe techniques used in our computer program 
for algebraically expanding Eq. I|19|l in orders of r. For 
a particular order of r, an expansion of the exponentials 
of U±(t) [Eq. I)17p] will result in the trace of a sum of 
products of H ± , and H B . The H± factors [Eq. (TSJ}] can 
all be distributed to yield a sum of traces of products of 
Ti. B , and H An = \ iZ n A n I nz . Each one of the Hamil- 
tonian pieces, Tt An , and Ti B , involves sums over nuclear 
site indices. We can decrease the number of independent 
indices by one (not a major improvement, but enough to 
push this method one step further) by noting that 



\Ti B , H An ~\ — — ^ b nm (A n — A m )I n _I„ 



(CI) 



which takes away the independence of H An, s index. Let 
TL C = [H B , TL 1 . We can convert our sum of traced 
products of Ji An and Tt B into a sum of traced products 
of H and Ti B and a single TL C with one less indepen- 
dent summation index. The rest of this paragraph will 
briefly explain how this is done. Consider all of the terms 
that involve a certain number of 7i s 's and 7i An, s. For 
any term, we can permute the Ti An, s and Tt B, s however 
we want at the expense of creating terms in which an 
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H B TtAn combination is replaced by TL C (which is the 
type of term we want anyway) . Consider the set of terms 
that all have a certain number of TCau factors and a cer- 
tain number of 7i B factors. As long as the coefficients 
of such a set add up to zero, we can make the permuta- 
tions cancel out, leaving only the terms with Ti c and one 
less independent summation index. By the symmetry in 
e ±iH t Q f gq g ^ an( j . these coefficients add up 
to zero for all such sets in any order of r. 

Consider pulling the summations over nuclear site in- 
dices, not only outside of each product, but also outside 
of the trace (since the trace of a sum is equal to the sum 
of the traces). The trace of these spin operators involving 
particular nuclear sites is straightforward to compute. In 
fact, all one really needs to know to perform the trace is 
which indices are the same and which ones are different 
(as well as the magnitude of the nuclear spins). Then one 
can compute the trace as follows. If O n is some operator 
(a product of spin operators I nzi I n +, In-) that operates 
on the spin of site n, then 

-Tt{O l O i ...O k )-— (C2) 

where M n is the number of states of the n th spin. To 
compute Tr (O n ) — J2i (*l On \i) we simply use 



I±\I,m) 
I z \I,m) 



a/J(I + 1) - m(m ± 1) | J", mil) ,(C3) 
m\I,m). (C4) 



Of course, in order to obtain an algebraic expression 
for a particular order of t, we do not want to consider 
the sum of all possible nuclear site indices. In fact, in 
order to obtain this expression, we should not even need 
to know how many nuclear spins there are. Instead, we 
simply consider all of the possible "index configurations." 
An "index configuration" indicates which indices are the 
same and which are different. As indicated above, this is 
enough information to allow one to evaluate the trace of 
the product of spin operators as a rational number [the 
square roots in Eq. (|C3(1 will always be repeated as factors 
an even number of times in the trace]. Additionally, we 
need only consider "index configurations" that do not 
obviously result in a zero trace for any of the indices. 
In particular, we only need to consider configurations in 
which there are the same number of i+ as J_ operators 
for each index, and we can ignore any configuration in 
which the only operator for any index is an odd power of 

Once the operators and traces are taken care of and 
we have exhausted all "index configurations," we are left 
with terms that have factors of A n 's and b nm 's, some 
with the same nuclear index labels, and some with differ- 
ent labels. These labels are to be summed over distinctly 
(i.e., different labels should never be assigned the same 
index). The labels are arbitrary, so when checking for 
"like" terms, in order to provide a compact expression 
for the result, we should consider that the labels can be 



permuted. It would be a waste of computation time to al- 
ways run through all permutations of labels when check- 
ing to see if two terms are like terms. Instead, the prob- 
lem is transformed into that of checking for graph iso- 
morphisms, which is a well-studied computational prob- 
lem. The program uses a code called Nauty written by 
Brendan McKay (http://cs.anu.edu.au/people/bdm) for 
checking for graph isomorphisms. 

The A n and b nm factors are transformed into a graph 
in the following way. Each distinct label is represented 
by a vertex in the graph. A factor of A\ is represented 
by coloring the vertex i to indicate the power p. A factor 
of bij is represented by an edge from vertex i to vertex 
j. If bij is raised to some power (other than 1) then it is 
represented by making an edge from both i and j to some 
special vertex that is colored in such a way as to indicate 
that it represents the appropriate power of bij (this is 
just a trick to effectively "color" edges). Then the Nauty 
program takes this graph and gives a canonical labeling 
for the vertices. If the program comes to another term 
of An and b nm factors that are equivalent to this one up 
to a permutation of labels, Nauty will give it the same 
canonical labeling so that the two terms can be combined. 



APPENDIX D: COMPUTATION OF 
TAU-EXPANSION COEFFICIENTS 

Given an algebraic expression for a particular order, in 
t, of ve(t), using the procedure described in Appendix 
IO we need to compute numerical coefficients for specific 
applications (with explicit values of A n and b nm ). For 
a problem that has N nontrivial nuclei and an algebraic 
expression that has m terms and k summation labels, 
this will in general take 0(mN k ) time to compute. This 
is a huge improvement over direct methods of evolving 
the state or diagonalizing the Hamiltonian since there 
are 0[M = (21 + 1) N ] states of the system. So this 
expansion already saves us from calculations that grow 
exponentially with N . However, when N is on the order 
of 1 000, which is actually a low estimate for considered 
applications, and we want the C(t 10 ) term which has 
k = 5, N k is still quite large (10 15 ). The situation can be 
drastically improved once more if we are willing to make 
another approximation which is very reasonable. We can 
set a threshold for values of |6„ m |, taking any \b nm \ be- 
low this threshold to be zero. Since \b nm \ generally de- 
creases as nuclei are further separated, this amounts to a 
near neighbor approximation (also exploited by our clus- 
ter method). We can then adjust this threshold until 
convergence is achieved. 

Let us say that our threshold is set such that, on aver- 
age, each nuclei couples to L other nuclei. If we assume 
that all of the terms of the expression to be calculated 
are connected, not disjoint, graphs when representing la- 
bels as vertices and factors of b nm as edges (as discussed 
in Appendix iDl in the context of checking for like terms), 
then this approximation reduces the computation time 
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to 0(mNL k ~ 1 ). For the first label, there are N possible 
nuclei to which it can be assigned, but each additional 
label must be one of the O(L) vertices connected by an 
edge to a previous vertex/label. 

Not all of the terms, in general, have these connected 
graph representations, however. Some of them are dis- 
joint. This does not create a real problem though. These 
"disjoint" terms can be factored and dealt with sepa- 
rately. As a simple example, 

distinct 

h v hki = E^E &fc;_2 E & y 

i,j,k,l i^j k^l ijtj 

distinct 

-4 ]T bijb jk , (Dl) 

i,j,k 

where the summation on the left (and one of them on 
the right) is over distinct indices as indicated above the 
summation symbol. The second term on the right has the 



factor of 2 because it includes bijbji 



b\, (since b r 



bmn)- The last term gets a factor of 2 from bjk — b^j and 
another factor of 2 from 



E 



bijbik 



E 



E b iO b jk, 



(D2) 



in which we relabeled indices. 

Our computer program, described in Sec. IV B II au- 
tomatically performs these conversions as a preprocess 
before being fed to the program that performs the nu- 
merical calculations. In addition to this conversion for 
"disjoint" terms, the preprocessor performs a simple uni- 
directionally cascading factorization, for example, of the 
form a[c + die + /)] but not (a + b)(c + d), which reduces 
the number of multiplications that must be performed in 
the calculation. This factorization does not change the 
time complexity of the calculation, which is 0(mNL k ^ 1 ), 
but it does provide a noticeable speedup. 



APPENDIX E: DIPOLAR PERTURBATION AND 
CLUSTER SIZE 

In this appendix we prove that the contribution of 
a cluster of size k is 0(X k ) with respect to the dipo- 
lar perturbation expansion of Sec. IV CI at least when 
t -C max(& nm ) _1 . A direct application of quantum per- 
turbation theory to the Halm echo yields Eq. I|53l) , giving 
a sum of oscillations of amplitude Cijki — Fiju / and 



frequency a;. 



,(°) 



XcoL 



Hjkl — ^ijkl T ^ijkl- 

First, let us only consider the nuclei involved in a term 
of Fijki (the numerator of Cijki ) ■ The way in which nuclei 
"get involved" in F^ki is via the sum over nuclear pairs in 
HI = {H B . This gives potentially two nuclei for one A, 
but we need to do better than that. Substituting \k±) on 
the right-hand side of Eq. I|47|l with the full expression 
will give successive orders of A. When one considers a 



particular term of u° \k±) with the above recursion for- 
mula in mind, one will go from state |Z°) to state 
with some number of intermediate Tt° eigenstates in be- 
tween where each state is connected to the next via H! . 
The second term in the numerator of Eq. i|47[l simply cor- 
responds to going from state |/°) to through some 
number of intermediate states and then from \ k } back to 
|fc°) through some more intermediate states. The point 
is simply that one goes from one state to a different state 
only via W [Eq. JTSJl]. With this in mind, Eq. JSSJl indi- 
cates that one must go from state |«°) to to to 
|/°) and back to state with some number of states in 
between where consecutive states are connected via H! . 
Because we must come back to state for every lower- 
ing operator we pass through, we must also pass through 
a raising operator for the same spin. We now see that we 
cannot actually get two nuclei for one A; for each nucleus 
involved in a given term of fykj, there is at least one 
order of A. 

The denominator, D+jki, works in a similar way, except 
that instead of one cycle from state back to state 
one has one of these going from to another 
going from to another going from |fc°) to 
and another going from |Z°) to But the result is the 
same; because these are all cyclic, there must be at least 
one order of A for each nucleus involved. One then needs 
to perform the Taylor expansion of D~^. Let g(X) — 

Dijki(X) and /(A) = ^jy- The Taylor series expansion is 



/(A) = £/ (n) (0) 



A" 



(El) 



n=0 



where /<» (A) may be obtained recursively with 



\(gW) k 



7 (m+l) 



(A) 



g ( m )(A) 
: (.9(A)) fe+1 



(E2) 



and then we note that [ff(0)] fc = 1 so there will not be 
anything in the denominators of Eq. i|E2|) to worry about 
in the context of Eq. i|Elfl . Examining the numerators 
in Eq. (|E2f) . we note that for each derivative of /(A), we 
will lose at most one order of A. But in Eq. I|E1|I for each 
derivative there is a factor of A provided to compensate. 
If we write the A expansion of 5(A) as 



7(A) 



5«A 



(E3) 



then, when we make the transformation to /(A) = 
l/g(X), each instance of g n will be accompanied by at 
least n orders of A. Therefore, after Taylor expanding 
the inverse of Dijki, we preserve the property that there 
are at least as many orders of A as the number of nuclei 
involved in a given term. 

We have dealt with the numerator and denominator of 
the amplitude, but we must still consider the frequency. 
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We can rewrite the exponential in Eq. I|53[) as 

e -imj k ir = e -^l lT exp (_ iAa ^. fe;T ) ( E4 ) 

= e IJ ~\ • ( E5 ) 

71=0 

Thus we have transformed the perturbative part of the 
frequency into an amplitude (with a time dependence as 
well) . We can now use a similar argument that was used 
for Fijki and Dijki- Each term of Eq. (|5*^)l will cycle from 
some state back to the same state with H' connecting 
each intermediate state. As a result, there must be at 
least as many orders of A in a given term of Xu>^ kl as there 
are nuclei involved. There will, however, be one extra 
&nm-type factor without an accompanying l/(A n — A m ). 
This gets combined with r and will be small as long as 
max (6„ m )r <C 1. All of this can be exponentiated in the 
series of Eq. l)E5jl , but that can only increase orders of A 
without increasing the number of nuclei involved. 

There is one final consideration. How many nuclei are 

involved in exp f — iw^rj ? From Eqs. (|58|l and l|49|) . we 

see that this will include the nuclei that differ between 
states and along with the nuclei that differ be- 



tween states |fc°) and |^°). Since a term of F^ki will 
require linking state to to |/c°) to |^°) and back 
to via H! (through some number of nuclear states), 
it must involve at least that many nuclei. So this factor 
will neither increase the number of involved nuclei nor 
change the number of A factors. 

We therefore conclude that, using an approximation 
that assumes max(6„ m )T 1, there will be at least as 
many orders of A in a term of the A-expansion of Eq. I|53|) 
as there are nuclei involved. Therefore a cluster contri- 
bution of size k will be of 0(X k ) or a higher order of 
A. This justifies the cluster expansion in the regime in 
which c nm 1 and max(6„ m )r <C 1. To be more accu- 
rate and to relax the max (6„ m )r <C 1 constraint, as long 
as max(6 nm )T < 1 we have shown that clusters of size k 
will give a contribution of order A fc_1 . Because it is con- 
sistent for different cluster sizes, k, and because clusters 
of size one do not contribute, this reduction by an order 
in A does not have an impact upon the validity of the 
cluster expansion in any of its various forms (i.e., prod- 
uct form or exponentiated form). For simplicity of the 
discussion, however, we will regard max (b nm )T simply as 
another order of A. 
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